High-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method and system based on moment method and equivalent principle method

Abstract

The application relates to a high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method and system based on a moment method and an equivalent principle method, high-order line basis functions are introduced into a rotationally symmetric body algorithm based on the equivalent principle method, adaptive high-order matrix filling, loop reduction and boundary division are completed based on a deep calculation unit architecture, and an LU matrix inversion operator under the architecture is proposed; first, a rotationally symmetric target is geometrically modeled and equivalent surface modeled, an EPA-BoR equation set is constructed, a high-order U matrix and a transfer matrix are derived; then, each matrix of the high-order EPA-BoR equation set is generated based on a heterogeneous parallel calculation strategy; after each MPI process generates corresponding scattering operators and translation operators, electromagnetic scattering solving is completed through inter-process transmission calculation; and finally, the radar scattering cross section of the rotationally symmetric target is calculated according to current solving results. The application has high calculation precision, low memory occupation and greatly improved solving efficiency.

Description

Technical Field

[0001] This invention belongs to the field of electromagnetic computing technology, specifically relating to a high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method and system based on the method of moments and the equivalence principle. Background Technology

[0002] When solving electromagnetic problems, speed and accuracy are generally of primary concern. High-frequency methods are computationally fast and memory-efficient, but lack accuracy. Among low-frequency methods, the finite element method (FEM) can generate sparse matrices, but it requires the introduction of absorbing boundary conditions, resulting in numerous unknowns and mesh truncation errors. The method of moments (MoM), based on integration, offers high accuracy but incurs high memory and time costs. The Fast Multilevel Submethod (FMM) addresses this by achieving efficient solutions to electromagnetic problems while maintaining high accuracy. For electrically large rotationally symmetric targets, the Body of Revolution - Method of MoM (BoR-MoM) algorithm, based on the MoM, reduces the three-dimensional electromagnetic scattering problem to a two-and-a-half-dimensional problem by leveraging rotational symmetry, significantly reducing the number of unknowns and achieving a solution speed far exceeding that of the FMM.

[0003] In recent years, parallel acceleration techniques for higher-order rotationally symmetric body moment methods and high-performance computing architectures have provided a core technical path for high-precision and high-efficiency solutions to electromagnetic scattering of electrically large rotationally symmetric targets. For example, some studies have addressed the core pain point of traditional BoR-MoM in solving electrically large targets by proposing a higher-order BoR-MoM scheme within a parallel framework. This scheme improves the accuracy of electromagnetic flow fitting on the target surface by introducing higher-order basis functions, while supporting distributed expansion across multiple GPUs, providing a mature heterogeneous parallel technical path for fast and high-precision solutions to electrically large rotationally symmetric targets. Chinese patent CN 103279589 B discloses a simulation method for electromagnetic scattering characteristics of rotationally symmetric bodies based on matrix nested compression. This method addresses the core bottleneck of high storage and computational complexity and rapidly increasing memory consumption with increasing electrical size in traditional BoR-MoM methods for solving electrically large targets. It constructs a binary tree grouping system for segmenting busbars, combines low-rank matrix decomposition with layer-by-layer nested compression to reduce the dimensionality of the impedance matrix, and introduces high-oscillation integration techniques to accelerate the calculation of high-mode Green's functions. This provides an efficient numerical simulation scheme for low-memory, fast-convergence solutions for electrically large rotationally symmetric targets. Chinese patent CN 118797908 B discloses a heterogeneous parallel fast direct solution method for accelerating the solution of complex electromagnetic problems. This method addresses the problems of high complexity and low utilization of heterogeneous computing resources in traditional methods of moments. It completes matrix LU decomposition and core computation acceleration through a CPU-GPU collaborative strategy, and combines an adaptive cross-approximation algorithm to achieve matrix compression and dimensionality reduction. This provides an innovative heterogeneous parallel technology path for efficient and high-precision direct solutions to complex electromagnetic problems.

[0004] However, the aforementioned high-order rotationally symmetric body moment method, matrix nested compression solution scheme, and heterogeneous parallel direct solution method still have significant shortcomings when dealing with large-scale dynamic rotationally symmetric cluster targets. While the high-order BoR-MoM scheme within a parallel framework improves the accuracy of electromagnetic flow fitting through higher-order basis functions and supports multi-GPU distributed expansion, it fails to incorporate the characteristics of the equivalence principle algorithm. When dealing with multiple clustered targets, the data communication overhead between sub-targets is significant, easily affecting parallel efficiency and making it difficult to effectively reduce the scale of unknowns in multiple clustered targets. The BoR solution method based on matrix nesting compression, although significantly reducing the memory and computational overhead of solving rotationally symmetric bodies through low-rank compression, does not introduce higher-order basis functions and distributed parallel expansion, making it unsuitable for the efficient solution requirements of large-scale rotationally symmetric body clusters. The fast direct solution method based on a CPU-GPU heterogeneous architecture achieves general parallel acceleration and direct solution for complex electromagnetic problems, but its parallel strategy is not customized for the mode orthogonality and dimensionality reduction characteristics of rotationally symmetric bodies. Furthermore, when calculating the inverse of the LU decomposition, it cannot fully utilize idle topological grid computing units for solving the synchronous LU inverse matrix, and its solution efficiency needs further improvement. Summary of the Invention

[0005] The purpose of this invention is to provide a high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method and system based on the method of moments and the equivalence principle, which has high computational accuracy, low memory usage, and high solution efficiency.

[0006] The technical solution to achieve the purpose of this invention is: a high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle, comprising the following steps:

[0007] Step 1: Perform geometric and equivalent surface modeling for the rotationally symmetric targets: For the cluster of targets to be analyzed, perform independent geometric modeling for each rotationally symmetric target; based on the Huygens equivalence principle, construct a closed equivalent surface outside each rotationally symmetric target; discretize the generatrices of the rotationally symmetric targets and the equivalent surface.

[0008] Step 2: Construct the EPA-BoR equation system and derive the higher-order U matrix and transmission matrix: Establish integral equations for the modeled rotationally symmetric target and equivalent surface; construct the corresponding EPA-BoR equation system based on the method of moments; derive and establish the higher-order U matrix and transmission matrix based on the formulas for higher-order basis functions and higher-order impedance matrices; substitute the formulas for the higher-order U matrix and transmission matrix into the EPA-BoR equation system to construct the higher-order EPA-BoR equation system.

[0009] Step 3: Generate matrices for the high-order EPA-BoR equations based on a heterogeneous parallel computing strategy: For the established high-order EPA-BoR equations, use an MPI process to divide them into blocks. Each MPI process independently calculates the corresponding scattering and translation operators based on the incident field information. For the high-order matrices contained in the scattering and translation operators, use a DCU architecture to perform heterogeneous filling calculations, cyclic reduction of the solution results from each thread, and boundary partitioning. For the calculated high-order impedance matrix, use a DCU architecture to perform LU decomposition and inversion calculations.

[0010] Step 4: After each MPI process generates its corresponding scattering and translation operators, the electromagnetic scattering solution is completed through inter-process transfer calculation: the preliminary equivalent surface electromagnetic flow solution is obtained through MPI transfer calculation; the error is reduced by using an iterative solution method to obtain the final local solution for each subdomain;

[0011] Step 5: Calculate the radar cross section of the rotationally symmetric target based on the current solution.

[0012] Further, step 1 includes:

[0013] Step 1.1: Model each rotationally symmetric target in the cluster of targets to be analyzed independently. A rotationally symmetric target is a spatial geometric model formed by rotating any continuous curve without repeating points around a straight line. The continuous curve is the generatrix, and the straight line is the axis of rotation.

[0014] Step 1.2: Based on Huygens' equivalence principle, construct a closed equivalent surface that is coaxial with the target and also rotationally symmetric outside each rotationally symmetric target. The equivalent surface is spaced at a set distance from the target surface.

[0015] Step 1.3: Discretize the generatrices of the rotationally symmetric target and the equivalent surface. After introducing higher-order line basis functions, set the discretization length to three-tenths of the wavelength corresponding to the incident wave frequency.

[0016] Step 1.4: Number the points and line segments on the generatrix, and record the vector coordinates of the normal vector of each line segment and the three-dimensional coordinates of the midpoint of the line segment.

[0017] Further, step 2 includes:

[0018] Step 2.1, Establishment of integral equation: When the incident wave irradiates the rotationally symmetric target, an induced electromagnetic current is generated on the target surface. The induced electromagnetic current generates a scattered field in the far region. Based on the relationship between the incident wave and the surface electromagnetic current, an integral equation of the electromagnetic field is established. The electromagnetic current on the target surface is divided into two directions: tangential and normal.

[0019] Step 2.2, Establishment of the EPA-BoR equation set: Based on the principle of the method of moments, the constructed electromagnetic field integral equation is transformed into a matrix equation form; the matrix equations are integrated to construct the EPA-BoR equation set describing the relationship between the equivalent surface scattered current and the incident current.

[0020] Step 2.3: Based on the formulas for higher-order basis functions and higher-order impedance matrices, derive and establish higher-order U matrices and transmission matrices: The U matrix is ​​a tridiagonal matrix, which is extended to a pentagonal matrix after applying second-order basis functions; the transmission matrix includes the transmission matrix between the equivalent surface and the internal rotationally symmetric target, and the transmission matrix between different equivalent surfaces. The transmission matrix is ​​composed of KL transfer operators. The KL operator form is extended to a higher order and filled into the original transmission matrix to form a higher-order transmission matrix.

[0021] Step 2.4: Substitute the higher-order impedance matrix and the derived higher-order U matrix and higher-order transmission matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations.

[0022] Furthermore, step 2.2 is detailed as follows:

[0023] calculate When there are rotationally symmetric targets =1,…, , =1,…, , , They are the first Scattered current density and magnetic flux density on the equivalent surface of a rotationally symmetric target. , The score is the first Incident current density and magnetic flux density on the equivalent surface of a rotationally symmetric target , They are the first Scattered current density and magnetic flux density on the equivalent surface of a rotationally symmetric target. It is the first Scattering operator for a rotationally symmetric target It is the first The rotationally symmetric target and the first The translation operator between the rotationally symmetric target equivalent surfaces, and the EPA-BoR equations are as follows:

[0024]

[0025] The matrix equation is as follows:

[0026]

[0027] in, For the first The current coefficient matrix of a rotationally symmetric target is the U matrix. For the first A rotationally symmetric K-operator describing the source field relationship. For the first An L-operator that describes the source field relationship using rotationally symmetric targets. For the first The impedance matrix of a rotationally symmetric target; It is the unit vector normal to the equivalent surface, pointing to the region where the field point is located.

[0028] Furthermore, step 2.3 is detailed as follows:

[0029] (1) Number the second-order basis functions. If the number of segments of the generatrix is ​​N, the second-order basis functions are numbered as the starting number 1. Then the first-order basis functions and the second-order basis functions are numbered alternately, and the second-order basis functions are numbered as the ending number 2N-1.

[0030] The expression for the second-order basis functions is as follows:

[0031]

[0032] in, , These represent the first and second segments of the piecewise basis function one, respectively. Indicates the second basis function; Represents the first in the parametric coordinate system part The tangent vector at the point; Indicates the coordinates of the parameters;

[0033] The expression for the elements of matrix U is as follows:

[0034]

[0035] in, To verify the index of the basis functions, represents the index of the basis functions in the current expansion. Indicates definition at position And the serial number is The test basis functions, Indicates definition at position And the serial number is The source basis functions;

[0036] Substituting the second-order basis function expression into the element expression of the U matrix, each element of the second-order U matrix is ​​the inner product of basis functions. When the two basis functions do not intersect and their inner product is zero, the U matrix is ​​expanded into a pentagonal matrix, as follows:

[0037]

[0038] in, The basis functions are represented sequentially as follows: 1, 2, 3, ..., 2N-3, 2N-2, 2N-1;

[0039] (2) Based on the different types of fields and flows, the KL transfer operators in the transfer matrix are divided into two categories:

[0040] electric field The operators relating to the interaction between current and magnetic current are respectively , ;

[0041] magnetic field The operators relating to the interaction between current and magnetic current are respectively and ;

[0042] in, , These represent the field point position vector and the source point position vector, respectively. Expressions and Consistent; Expressions and In comparison, free space conductivity Replace with free space permeability ;

[0043] Based on the method of moments, the matrix form of the KL transfer operator consists of four parts, which are the result of the interaction between the current in the tangential direction and the current in the circumferential direction of the busbar. The expressions of the sub-matrices are as follows:

[0044]

[0045]

[0046] in, Represents electric field With magnetic current The interaction matrix between them Represents electric field With current The interaction matrix between them Represents the first in the matrix OK Column elements; Indicates definition at position And the serial number is The source basis functions; Indicates definition at position And the serial number is The source basis functions; , These represent the surface element vector at the field point and the surface element vector at the source point, respectively. For free space wavenumber, Permeability in free space Let be the Euclidean linear distance between the field point and the source point. For the time-harmonic field scalar Green's function in three-dimensional free space, Angular frequency, The phase constant, Represents the imaginary unit; Represents the gradient operator; , These represent the surface divergence operators acting at the field point and the source point, respectively.

[0047] Substituting the expression for the second-order basis function into the transfer matrix, the matrix... OK The elements at each column are summed from the calculation results of each subdivision segment. Based on the properties of higher-order line basis functions, only the test basis function is used. and source basis functions The subdivision segment contains non-zero values; since the basis functions span two adjacent subdivision segments, in actual solution, the first... OK The elements in the column are the result of the interaction of four elements from two adjacent subdivision segments.

[0048] Further, step 3 includes:

[0049] Step 3.1: For the established high-order EPA-BoR equation set, use the MPI process to divide it into blocks: each MPI process independently reads the geometric and material parameters of the assigned BoR target and classifies the target as a metallic target or a dielectric target, and establishes different MPI task divisions; each MPI process independently calculates the corresponding scattering operator and translation operator based on the incident field information.

[0050] Step 3.2: For the higher-order matrices contained in the scattering and translation operators, a DCU architecture is used for heterogeneous filling calculations: No inter-process communication is required during the matrix filling stage; each process independently completes the filling of its assigned submatrix. Matrix filling is achieved by iteratively traversing the partitioned segments corresponding to the basis functions and test functions, as detailed below:

[0051] For the allocation Each test function and For a single MPI process with basis functions, the matrix elements corresponding to the computation are as follows:

[0052]

[0053] in, The serial number is the Fourier model number. , This represents the highest order of the Fourier module; Indicates the first The first Fourier mode OK Column matrix elements; Indicates the integral kernel The applied linear operator; This indicates the node number of the Gaussian integral of the test function. This indicates the node number of the Gaussian integral of the basis functions. Indicates the index of the Gaussian integration node of the integration kernel; This represents the total number of nodes in the Gaussian integral of the test function. This represents the total number of nodes in the Gaussian integral of the basis functions. This represents the total number of Gaussian integration nodes in the integration kernel; Indicates the first Gaussian integration weight coefficients corresponding to the Gaussian integration nodes of the test function. Indicates the first Gaussian integration weight coefficients corresponding to the Gaussian integration nodes of the basis functions. Indicates the first Gaussian integration weight coefficients corresponding to each Gaussian integration node of the integration kernel;

[0054] First, all basis functions required for matrix filling are generated and transferred to the accelerator's global memory. Then, the accelerator maps the integration calculation task to various threads based on the matrix dimensions and integration node information, where the grid dimensions are determined by... Configure the thread block dimension by Configure;

[0055] Step 3.3: Perform heterogeneous cyclic reduction and boundary partitioning on the solution results of each thread in the DCU: After each thread completes the calculation of its own integral value, iterative updates of the cyclic reduction are performed using a cyclic reduction scheme; after the reduction is completed, the boundary is re-partitioned based on the properties of higher-order basis functions, as follows:

[0056] After each thread calculates the integral value, the result is reduced to form the complete matrix elements. A cyclic reduction scheme is adopted, and a one-dimensional index is reassigned to the integral value of each thread. The iterative update of the cyclic reduction is completed according to the following formula:

[0057]

[0058] Among them, the one-dimensional index in the specification process The three dimensions of a thread block are as follows: ; This refers to the sequence number of the protocol round. Indicates the first The index in the round specification is The integral value, Indicates the first The index in the round specification is The integral value, Indicates the first The index in the round specification is The integral value;

[0059] Through iterative reduction by pairwise grouping round after round, the final result is achieved... After round reduction, the index is obtained. The integral value of 1 This is the final summation of the integral values ​​of each thread corresponding to the matrix element;

[0060] Based on the properties of second-order basis functions, when filling matrix elements at the end of reduction accumulation, to ensure that each DCU thread operates completely independently, each MPI process is responsible for calculating the matrix element values ​​on a segment of the busbar; the basis functions of adjacent processes are second-order. They will not intersect; the first segment of piecewise basis function one. The second segment of segmented basis function one, allocated to the previous process. The matrix is ​​assigned to the next process; after each process has calculated its corresponding block matrix, the matrix boundaries of two adjacent processes are overlapped and added together to obtain the result of the corresponding part of the original matrix.

[0061] Step 3.4: For the calculated high-order impedance matrix, perform LU decomposition and inversion calculation using a DCU architecture: In the direct solver based on LU decomposition, the three core computationally intensive operations, LU decomposition, triangular matrix inversion, and matrix multiplication, are all executed at the DCU end; the heterogeneous computing interface basic linear algebra subroutine library hipBLAS is used to accelerate LU decomposition and matrix multiplication.

[0062] For the triangular matrix inversion operation in the second step of parallel decomposition, an acceleration operator is designed to speed up the triangular matrix inversion part on the DCU, as follows:

[0063] Based on a uniform column-strip partitioning strategy with fixed block granularity, two types of triangular matrices are solved using delayed-correction recursive frameworks from right to left and from left to right, respectively; the triangular matrices to be inverted are uniformly denoted as... , The unit lower triangular matrix corresponding to LU decomposition , The non-singular upper triangular matrix corresponding to LU decomposition According to fixed block particle size Will 1-th order matrix Evenly divided into Each column of blocks, the inverse matrix to be solved. Adopted and Completely consistent column-bar block partitioning rules and initial values ​​set to identity matrix of order The unified delay correction recursive formula is as follows:

[0064]

[0065] In the formula Number the columns and blocks. Inverse matrix The Each column of bars, For matrix The Middle Each column block corresponds to The inverse matrix of the main diagonal sub-block. for An identity matrix of order 1. For the first Excluding the main diagonal pieces from each column of blocks. The strictly triangular submatrix after that, This is the set of inverse sub-blocks that have been solved in the recursive direction;

[0066] when When it is a unit lower triangular matrix, The set of column numbers to the right of the current column block. The recursive process is numbered by column block. from Execute from right to left until step 1; when When it is a non-singular upper triangular matrix, The set of column numbers to the left of the current column block. The recursive process is numbered by column block. From 1 to Execute from left to right;

[0067] At the thread mapping level, a two-dimensional Cartesian topology is used to construct the thread grid and thread blocks: first, construct... Two-dimensional thread block, The number of threads is allocated to the corresponding side length of the thread block, and then a corresponding two-dimensional grid dimension is constructed based on the number of matrix blocks, so that the two-dimensional coordinates of the threads and the row and column indices of the matrix form a one-to-one mapping relationship; in terms of block partitioning strategy design, a two-level block partitioning scheme combining coarse-grained and fine-grained approaches is adopted: the first level is... The coarse-grained columnar block partitioning decomposes the large-scale global matrix into sub-block tasks adapted to the shared memory capacity on the DCU chip; the second level is... Fine-grained micro-block partitioning, dividing each The sub-blocks are further broken down into The micro-computing units are mapped to a single thread for completion, and all recursive calculations within the micro-block are executed in the thread's private registers, adapting to the execution pipeline characteristics of the DCU.

[0068] Further, step 4 includes:

[0069] Step 4.1: Solve the initial equivalent surface electromagnetic flow solution through MPI transfer calculation: Each MPI process solves the initial equivalent surface current on the corresponding equivalent surface ES, which will be used as the initial solution of the iterative solver; the target initialization is completed in each MPI process to preserve the locality of electromagnetic interaction;

[0070] Step 4.2: Use iterative solution method to reduce error and find the final local solution of each subdomain: In the iterative solution stage, the generalized minimum residual (GMRES) method is adopted to update the local solution of each subdomain through scattering operator and translation operator until the relative residual norm converges to below the preset threshold.

[0071] Further, step 5 includes:

[0072] Step 5.1: After completing the iterative solution, set the equivalent surface electromagnetic currents obtained by convergence in each subdomain as independent Huygens equivalent radiation sources.

[0073] Step 5.2: For each equivalent radiation source, use the far-field radiation integral formula to calculate the far-field scattered electric field vector of the equivalent radiation source in the specified observation direction; then, superimpose the far-field scattered electric field vectors generated by all equivalent radiation sources in the same observation direction to obtain the total scattered electric field vector of the entire rotationally symmetric target in that observation direction.

[0074] Step 5.3: Based on the standard definition of radar cross section, calculate the radar cross section value of the target in the observation direction, defined as the ratio of the power density scattered by the target at the receiver per unit solid angle to the incident power density of the incident wave on the target surface. times.

[0075] A high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation system based on the method of moments and the equivalence principle is provided. This system is used to implement the aforementioned high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle. The system includes:

[0076] Geometric modeling and equivalent surface modeling module for rotationally symmetric targets: For the cluster of targets to be analyzed, each rotationally symmetric target is independently geometrically modeled; based on Huygens' equivalence principle, a closed equivalent surface is constructed outside each rotationally symmetric target; the generatrices of the rotationally symmetric targets and the equivalent surface are discretized.

[0077] The module for constructing higher-order EPA-BoR equations establishes integral equations for the modeled rotationally symmetric target and equivalent surface; constructs the corresponding EPA-BoR equations based on the method of moments; derives and establishes the higher-order U matrix and transmission matrix based on the formulas for higher-order basis functions and higher-order impedance matrices; and substitutes the formulas for the higher-order U matrix and transmission matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations.

[0078] The matrix generation modules for the higher-order EPA-BoR equations are as follows: For the established higher-order EPA-BoR equations, the MPI process is used to divide them into blocks. Each MPI process independently calculates the corresponding scattering operator and translation operator based on the incident field information. For the higher-order matrices contained in the scattering operator and translation operator, the DCU architecture is used to perform heterogeneous filling calculations, cyclic reduction of the solution results of each thread, and boundary partitioning. For the calculated higher-order impedance matrix, the DCU architecture is used to perform LU decomposition and inversion calculations.

[0079] The electromagnetic scattering solution module generates its own corresponding scattering and translation operators in each MPI process, and then completes the electromagnetic scattering solution through inter-process transfer calculation: the preliminary equivalent surface electromagnetic flow solution is obtained through MPI transfer calculation; the error is reduced by using an iterative solution method, and the final local solution of each subdomain is obtained;

[0080] The radar cross section calculation module calculates the radar cross section of a rotationally symmetric target based on the current solution results.

[0081] A mobile terminal includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle.

[0082] Compared with the prior art, the significant advantages of this invention are:

[0083] (1) A new architecture for inverting the LU matrix is ​​proposed. First, the rotationally symmetric target is geometrically modeled and its equivalent surface modeled. Then, the EPA-BoR equations are constructed to derive the higher-order U matrix and the transfer matrix. Then, the matrix generation of the higher-order EPA-BoR equations is completed based on the heterogeneous parallel computing strategy. After each MPI process generates its corresponding scattering operator and translation operator, the electromagnetic scattering solution is completed through inter-process transfer calculation. The calculation accuracy is high, the memory usage is low, and the solution efficiency is greatly improved.

[0084] (2) The higher-order line basis functions are introduced into the rotational symmetric body algorithm EPA-BoR based on the equivalence principle method, and the adapted higher-order matrix filling, cyclic reduction and boundary partitioning are completed based on the deep computing unit (DCU) architecture. It has good parallel scalability and computational efficiency, and is particularly suitable for heterogeneous accelerated computing platforms.

[0085] (3) Introducing higher-order line basis functions effectively reduces the matrix size and memory usage while ensuring the solution accuracy. Based on the DCU architecture, the adaptation of higher-order matrix filling, cyclic reduction and boundary partitioning is completed. At the same time, a LU matrix inversion operator under the architecture is proposed, which significantly reduces memory requirements and communication overhead. Attached Figure Description

[0086] Figure 1 This is a flowchart of the high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle of the present invention.

[0087] Figure 2 This is a schematic diagram showing the dimensions and positions of the double cones in an embodiment of the present invention.

[0088] Figure 3 This is a schematic diagram showing the distribution and specific numbering of the second-order linear basis functions on the busbar in this invention.

[0089] Figure 4 This is a schematic diagram illustrating the specific iterative process of LU matrix inversion calculation based on the DCU architecture of this invention.

[0090] Figure 5 This is a schematic diagram comparing the calculated radar cross section (RCS) of the double cone with the calculation results of the simulation software FEKO in an embodiment of the present invention.

[0091] Figure 6 This is a schematic diagram comparing the unknowns of the RCS calculation for five conical targets at different frequencies using the second-order EPA-BoR algorithm and the first-order algorithm in an embodiment of the present invention.

[0092] Figure 7 This is a schematic diagram comparing the measured time consumption of the second-order EPA-BoR algorithm and the first-order algorithm for calculating RCS of five conical targets at different frequencies in an embodiment of the present invention. Detailed Implementation

[0093] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0094] This invention proposes a heterogeneous parallel electromagnetic simulation method for high-order rotationally symmetric targets based on the method of moments (MoM) and the equivalence principle (EPP). It introduces high-order line basis functions into the rotationally symmetric body algorithm (EPA-BoR) based on the EEPP, and utilizes a deep computational unit (DCU) architecture to perform adapted high-order matrix filling, cyclic reduction, and boundary partitioning. A LU matrix inversion operator under this architecture is also proposed. The method first performs geometric and equivalent surface modeling of the rotationally symmetric target; secondly, it constructs the EPA-BoR equations, deriving the high-order U matrix and transmission matrix; then, it generates the matrices of the high-order EPA-BoR equations using a heterogeneous parallel computing strategy; after each MPI process generates its corresponding scattering and translation operators, electromagnetic scattering is solved through inter-process transfer calculations; finally, the radar cross-section of the rotationally symmetric target is calculated based on the current solution. This invention offers high computational accuracy, low memory usage, and significantly improved solution efficiency.

[0095] Combination Figure 1 This invention presents a high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle. It introduces high-order line basis functions into the Equivalence Principle Algorithm - Body of Revolution (EPA-BoR) (EPA-BoR), leveraging the superior properties of high-order basis functions to achieve efficient analysis of the electromagnetic scattering characteristics of rotationally symmetric targets, significantly improving solution efficiency. Furthermore, based on a Deep Computing Unit (DCU) architecture, it completes adapted high-order matrix filling, cyclic reduction, and boundary partitioning, and proposes an LU matrix inversion operator under this architecture. The specific steps include:

[0096] Step 1: Perform geometric modeling and equivalent surface modeling for the rotationally symmetric target: such as... Figure 2 As shown, the model consists of two conical targets, each with a base radius of 1m, a height of 2m, and a distance of 8m. Each target element is geometrically modeled independently. Based on the Huygens equivalence principle, a closed equivalent surface is constructed outside each rotationally symmetric target. The generatrices of the rotationally symmetric body and the equivalent surface are discretized, as detailed below:

[0097] Step 1.1: Model each rotationally symmetric target in the cluster of targets to be analyzed independently. A rotationally symmetric target is a spatial geometric model formed by rotating any continuous curve without repeating points around a straight line. The continuous curve is the generatrix, and the straight line is the axis of rotation.

[0098] Step 1.2: Based on Huygens' equivalence principle, construct a closed equivalent surface that is coaxial with the target and also rotationally symmetric outside each rotationally symmetric target. The equivalent surface is spaced at a set distance from the target surface.

[0099] Step 1.3: Discretize the generatrices of the rotationally symmetric target and the equivalent surface. After introducing higher-order line basis functions, set the discretization length to three-tenths of the wavelength corresponding to the incident wave frequency.

[0100] Step 1.4: Number the points and line segments on the generatrix, and record the vector coordinates of the normal vector of each line segment and the three-dimensional coordinates of the midpoint of the line segment.

[0101] Step 2: Construct the EPA-BoR equations and derive the higher-order U matrix and transfer matrix: Establish integral equations for the modeled rotationally symmetric target and equivalent surface; construct the corresponding EPA-BoR equations based on the method of moments; derive and establish the higher-order U matrix and transfer matrix based on the formulas for higher-order basis functions and higher-order impedance matrices; substitute the formulas for the higher-order U matrix and transfer matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations, as follows:

[0102] Step 2.1, Establishment of integral equation: When the incident wave irradiates the rotationally symmetric target, an induced electromagnetic current is generated on the target surface. The induced electromagnetic current generates a scattered field in the far region. Based on the relationship between the incident wave and the surface electromagnetic current, an integral equation of the electromagnetic field is established. The electromagnetic current on the target surface is divided into two directions: tangential and normal.

[0103] Step 2.2, Establishment of the EPA-BoR Equations: Based on the principle of the method of moments, the constructed electromagnetic field integral equations are transformed into matrix equations; the matrix equations are then integrated to construct the EPA-BoR equations describing the relationship between the equivalent surface scattered current and the incident current, as follows:

[0104] calculate When there are rotationally symmetric targets =1,…, , =1,…, , , They are the first Scattered current density and magnetic flux density on the equivalent surface of a rotationally symmetric target. , The score is the first Incident current density and magnetic flux density on the equivalent surface of a rotationally symmetric target , They are the first Scattered current density and magnetic flux density on the equivalent surface of a rotationally symmetric target. It is the first Scattering operator for a rotationally symmetric target It is the first The rotationally symmetric target and the first The translation operator between the rotationally symmetric target equivalent surfaces, and the EPA-BoR equations are as follows:

[0105]

[0106] The matrix equation is as follows:

[0107]

[0108] in, For the first The current coefficient matrix of a rotationally symmetric target is the U matrix. For the first A rotationally symmetric K-operator describing the source field relationship. For the first An L-operator that describes the source field relationship using rotationally symmetric targets. For the first The impedance matrix of a rotationally symmetric target; It is the unit vector normal to the equivalent surface, pointing to the region where the field point is located.

[0109] Step 2.3: Based on the formulas for higher-order basis functions and higher-order impedance matrices, derive and establish the higher-order U matrix and the transmission matrix: The U matrix is ​​a tridiagonal matrix, which is extended to a pentagonal matrix after applying second-order basis functions; the transmission matrix includes the transmission matrix between the equivalent surface and the internal rotationally symmetric target, and the transmission matrix between different equivalent surfaces. The transmission matrix is ​​composed of the KL transfer operator. The KL operator form is extended to a higher order and filled into the original transmission matrix to form the higher-order transmission matrix, as detailed below:

[0110] (1) such as Figure 3 As shown, the second-order basis functions are numbered. If the number of segments of the busbar is N, the second-order basis functions are numbered as the starting number 1, and then the first-order basis functions and the second-order basis functions are numbered alternately, with the second-order basis functions as the ending number 2N-1.

[0111] The expression for the second-order basis functions is as follows:

[0112]

[0113] in, , These represent the first and second segments of the piecewise basis function one, respectively. Indicates the second basis function; Represents the first in the parametric coordinate system part The tangent vector at the point; Indicates the coordinates of the parameters;

[0114] The expression for the elements of matrix U is as follows:

[0115]

[0116] in, To verify the index of the basis functions, represents the index of the basis functions in the current expansion. Indicates definition at position And the serial number is The test basis functions, Indicates definition at position And the serial number is The source basis functions;

[0117] Substituting the second-order basis function expression into the element expression of the U matrix, each element of the second-order U matrix is ​​the inner product of basis functions. When the two basis functions do not intersect and their inner product is zero, the U matrix is ​​expanded into a pentagonal matrix, as follows:

[0118]

[0119] in, The basis functions are represented sequentially as follows: 1, 2, 3, ..., 2N-3, 2N-2, 2N-1;

[0120] (2) Based on the different types of fields and flows, the KL transfer operators in the transfer matrix are divided into two categories:

[0121] electric field The operators relating to the interaction between current and magnetic current are respectively , ;

[0122] magnetic field The operators relating to the interaction between current and magnetic current are respectively and ;

[0123] in, , These represent the field point position vector and the source point position vector, respectively. Expressions and Consistent; Expressions and In comparison, free space conductivity Replace with free space permeability ;

[0124] Based on the method of moments, the matrix form of the KL transfer operator consists of four parts, which are the result of the interaction between the current in the tangential direction and the current in the circumferential direction of the busbar. The expressions of the sub-matrices are as follows:

[0125]

[0126]

[0127] in, Represents electric field With magnetic current The interaction matrix between them Represents electric field With current The interaction matrix between them Represents the first in the matrix OK Column elements; Indicates definition at position And the serial number is The source basis functions; Indicates definition at position And the serial number is The source basis functions; , These represent the surface element vector at the field point and the surface element vector at the source point, respectively. For free space wavenumber, Permeability in free space Let be the Euclidean linear distance between the field point and the source point. For the time-harmonic field scalar Green's function in three-dimensional free space, Angular frequency, The phase constant, Represents the imaginary unit; Represents the gradient operator; , These represent the surface divergence operators acting at the field point and the source point, respectively.

[0128] Substituting the expression for the second-order basis function into the transfer matrix, the matrix... OK The elements at each column are summed from the calculation results of each subdivision segment. Based on the properties of higher-order line basis functions, only the test basis function is used. and source basis functions The subdivision segment contains non-zero values; since the basis functions span two adjacent subdivision segments, in actual solution, the first... OK The elements in the column are the result of the interaction of four elements from two adjacent subdivision segments.

[0129] Step 2.4: Substitute the higher-order impedance matrix and the derived higher-order U matrix and higher-order transmission matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations.

[0130] Step 3: Generate matrices for the high-order EPA-BoR equations based on a heterogeneous parallel computing strategy: For the established high-order EPA-BoR equations, use MPI processes to divide them into blocks. Each MPI process independently calculates the corresponding scattering and translation operators based on the incident field information. For the high-order matrices contained in the scattering and translation operators, use a DCU architecture for heterogeneous filling calculations, loop reduction of the solution results from each thread, and boundary partitioning. For the calculated high-order impedance matrix, use a DCU architecture for LU decomposition and inversion calculations, as detailed below:

[0131] Step 3.1: For the established high-order EPA-BoR equation set, use the MPI process to divide it into blocks: each MPI process independently reads the geometric and material parameters of the BoR target assigned to it, and classifies the target as a metallic target or a dielectric target, thereby establishing different MPI task divisions; each MPI process independently calculates the corresponding scattering operator and translation operator based on the incident field information.

[0132] Step 3.2: For the higher-order matrices contained in the scattering and translation operators, a DCU architecture is used for heterogeneous filling calculations: No inter-process communication is required during the matrix filling stage; each process independently completes the filling of its assigned submatrix. Matrix filling is achieved by iteratively traversing the partitioned segments corresponding to the basis functions and test functions, as detailed below:

[0133] For the allocation Each test function and For a single MPI process with basis functions, the matrix elements corresponding to the computation are as follows:

[0134]

[0135] in, The serial number is the Fourier model number. , This represents the highest order of the Fourier module; Indicates the first The first Fourier mode OK Column matrix elements; Indicates the integral kernel The applied linear operator; This indicates the node number of the Gaussian integral of the test function. This indicates the node number of the Gaussian integral of the basis functions. Indicates the index of the Gaussian integration node of the integration kernel; This represents the total number of nodes in the Gaussian integral of the test function. This represents the total number of nodes in the Gaussian integral of the basis functions. This represents the total number of Gaussian integration nodes in the integration kernel; Indicates the first Gaussian integration weight coefficients corresponding to the Gaussian integration nodes of the test function. Indicates the first Gaussian integration weight coefficients corresponding to the Gaussian integration nodes of the basis functions. Indicates the first Gaussian integration weight coefficients corresponding to each Gaussian integration node of the integration kernel;

[0136] First, all basis functions required for matrix filling are generated and transferred to the accelerator's global memory. Then, the accelerator maps the integration calculation task to various threads based on the matrix dimensions and integration node information, where the grid dimensions are determined by... Configure the thread block dimension by Configure;

[0137] Step 3.3: Perform heterogeneous cyclic reduction and boundary partitioning on the solution results of each thread in the DCU: After each thread completes the calculation of its own integral value, iterative updates of the cyclic reduction are performed using a cyclic reduction scheme; after the reduction is completed, the boundary is re-partitioned based on the properties of higher-order basis functions, as follows:

[0138] After each thread calculates its integral value, the result must be reduced to form the complete matrix elements. Since a large number of integral values ​​are accumulated simultaneously, direct accumulation can lead to severe read / write conflicts. Therefore, a cyclic reduction scheme is proposed, which reassigns a one-dimensional index to the integral value of each thread and performs iterative updates of the cyclic reduction according to the following formula:

[0139]

[0140] Among them, the one-dimensional index in the specification process The three dimensions of a thread block are as follows: ; This refers to the sequence number of the protocol round. Indicates the first The index in the round specification is The integral value, Indicates the first The index in the round specification is The integral value, Indicates the first The index in the round specification is The integral value;

[0141] Through iterative reduction by pairwise grouping round after round, the final result is achieved... After round reduction, the index is obtained. The integral value of 1 This is the final summation of the integral values ​​of each thread corresponding to the matrix elements; after adopting this method, the total number of reduction rounds is reduced from... Down to ;

[0142] Based on the properties of second-order basis functions, when filling matrix elements at the end of reduction accumulation, to ensure that each DCU thread operates completely independently, each MPI process is responsible for calculating the matrix element values ​​on a segment of the busbar; the basis functions of adjacent processes are second-order. They will not intersect; the first segment of piecewise basis function one. The second segment of segmented basis function one, allocated to the previous process. The matrix is ​​assigned to the next process; after each process has calculated its corresponding block matrix, the matrix boundaries of two adjacent processes are overlapped and added together to obtain the result of the corresponding part of the original matrix.

[0143] Step 3.4: For the calculated high-order impedance matrix, LU decomposition and inversion are performed using a DCU architecture. In the LU decomposition-based direct solver, the three core computationally intensive operations—LU decomposition, triangular matrix inversion, and matrix multiplication—are all executed on the DCU. The Portable Heterogeneous Computing Interface Basic Linear Algebra Subroutine Library (hipBLAS) is used to accelerate LU decomposition and matrix multiplication. For the triangular matrix inversion operation in the second step of parallel decomposition, a dedicated operator is designed to accelerate the triangular matrix inversion part on the DCU, as follows:

[0144] Based on a uniform column-strip partitioning strategy with fixed block granularity, two types of triangular matrices are solved using delayed-correction recursive frameworks from right to left and from left to right, respectively; the triangular matrices to be inverted are uniformly denoted as... , The unit lower triangular matrix corresponding to LU decomposition , The non-singular upper triangular matrix corresponding to LU decomposition According to fixed block particle size Will 1-th order matrix Evenly divided into Each column of blocks, the inverse matrix to be solved. Adopted and Completely consistent column-bar block partitioning rules and initial values ​​set to identity matrix of order The unified delay correction recursive formula is as follows:

[0145]

[0146] In the formula Number the columns and blocks. Inverse matrix The Each column of bars, For matrix The Middle Each column block corresponds to The inverse matrix of the main diagonal sub-block. for An identity matrix of order 1. For the first Excluding the main diagonal pieces from each column of blocks. The strictly triangular submatrix after that, This is the set of inverse sub-blocks that have been solved in the recursive direction;

[0147] like Figure 4 As shown, Procsx and Procsy are the first and second dimensions of the DCU's two-dimensional Cartesian topological mesh, respectively. When it is a unit lower triangular matrix, The set of column numbers to the right of the current column block. The recursive process is numbered by column block. from Execute from right to left until step 1; when When it is a non-singular upper triangular matrix, The set of column numbers to the left of the current column block. The recursive process is numbered by column block. From 1 to Execute from left to right;

[0148] At the thread mapping level, a two-dimensional Cartesian topology is used to construct the thread grid and thread blocks: first, construct... Two-dimensional thread block, The number of threads is allocated to the corresponding side length of the thread block, and then a corresponding two-dimensional grid dimension is constructed based on the number of matrix blocks, so that the two-dimensional coordinates of the threads and the row and column indices of the matrix form a one-to-one mapping relationship; in terms of block partitioning strategy design, a two-level block partitioning scheme combining coarse-grained and fine-grained approaches is adopted: the first level is... The coarse-grained columnar block partitioning decomposes the large-scale global matrix into sub-block tasks adapted to the shared memory capacity on the DCU chip; the second level is... Fine-grained micro-block partitioning, dividing each The sub-blocks are further broken down into The micro-computing units are mapped to a single thread for completion, and all recursive calculations within the micro-block are executed in the thread's private registers, adapting to the execution pipeline characteristics of the DCU.

[0149] This operator can compute the inverses of the L and U matrices simultaneously and in parallel. During the computation process, it can make full use of all threads in the thread grid constructed by the two-dimensional Cartesian topology, and keep a small number of threads idle during each iteration.

[0150] Step 4: After each MPI process generates its corresponding scattering and translation operators, the electromagnetic scattering solution is completed through inter-process transfer calculations: the preliminary equivalent surface electromagnetic flow solution is obtained through MPI transfer calculations; the error is reduced using an iterative solution method to obtain the final local solution for each subdomain, as detailed below:

[0151] Step 4.1: Solve the initial equivalent surface electromagnetic flow solution through MPI transfer calculation: Each MPI process solves the initial equivalent surface current on the corresponding equivalent surface ES, which will be used as the initial solution of the iterative solver; the target initialization is completed in each MPI process to preserve the locality of electromagnetic interaction;

[0152] Step 4.2: Use iterative solution method to reduce error and find the final local solution of each subdomain: In the iterative solution stage, the generalized minimum residual (GMRES) method is adopted to update the local solution of each subdomain through scattering operator and translation operator until the relative residual norm converges to below the preset threshold.

[0153] Step 5: Calculate the radar cross-section of the rotationally symmetric target based on the current solution results, as follows:

[0154] Step 5.1: After completing the iterative solution, set the equivalent surface electromagnetic currents obtained by convergence in each subdomain as independent Huygens equivalent radiation sources.

[0155] Step 5.2: For each equivalent radiation source, use the far-field radiation integral formula to calculate the far-field scattered electric field vector of the equivalent radiation source in the specified observation direction; then, superimpose the far-field scattered electric field vectors generated by all equivalent radiation sources in the same observation direction to obtain the total scattered electric field vector of the entire rotationally symmetric target in that observation direction.

[0156] Step 5.3: Based on the standard definition of radar cross section, calculate the radar cross section value of the target in the observation direction, defined as the ratio of the power density scattered by the target at the receiver per unit solid angle to the incident power density of the incident wave on the target surface. times.

[0157] This invention also provides a high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation system based on the method of moments and the equivalence principle. This system is used to implement the aforementioned high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle. The system includes:

[0158] Geometric modeling and equivalent surface modeling module for rotationally symmetric targets: For the cluster of targets to be analyzed, each rotationally symmetric target is independently geometrically modeled; based on Huygens' equivalence principle, a closed equivalent surface is constructed outside each rotationally symmetric target; the generatrices of the rotationally symmetric targets and the equivalent surface are discretized.

[0159] The module for constructing higher-order EPA-BoR equations establishes integral equations for the modeled rotationally symmetric target and equivalent surface; constructs the corresponding EPA-BoR equations based on the method of moments; derives and establishes the higher-order U matrix and transmission matrix based on the formulas for higher-order basis functions and higher-order impedance matrices; and substitutes the formulas for the higher-order U matrix and transmission matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations.

[0160] The matrix generation modules for the higher-order EPA-BoR equations are as follows: For the established higher-order EPA-BoR equations, the MPI process is used to divide them into blocks. Each MPI process independently calculates the corresponding scattering operator and translation operator based on the incident field information. For the higher-order matrices contained in the scattering operator and translation operator, the DCU architecture is used to perform heterogeneous filling calculations, cyclic reduction of the solution results of each thread, and boundary partitioning. For the calculated higher-order impedance matrix, the DCU architecture is used to perform LU decomposition and inversion calculations.

[0161] The electromagnetic scattering solution module generates its own corresponding scattering and translation operators in each MPI process, and then completes the electromagnetic scattering solution through inter-process transfer calculation: the preliminary equivalent surface electromagnetic flow solution is obtained through MPI transfer calculation; the error is reduced by using an iterative solution method, and the final local solution of each subdomain is obtained;

[0162] The radar cross section calculation module calculates the radar cross section of a rotationally symmetric target based on the current solution results.

[0163] The present invention also provides a mobile terminal, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle.

[0164] Example

[0165] For the analysis of electromagnetic scattering characteristics of multi-rotationally symmetric targets, traditional serial algorithms based on first-order basis functions face computational bottlenecks and storage limitations, making it difficult to meet the requirements of computational efficiency and accuracy in practical engineering.

[0166] This embodiment uses a double cone as an example, and the radar simulation parameters are set as follows: carrier frequency is 1 GHz, plane wave incident direction is perpendicular, and polarization mode is VV polarization. Figure 5As shown, the comparison between the calculated radar cross section (RCS) of the double cone and the calculation results from the simulation software FEKO verifies the accuracy of the calculation results of the parallel algorithm for high-order rotationally symmetric bodies based on the method of moments and the equivalence principle constructed in this invention. To further evaluate the parallel efficiency of the algorithm, this paper conducts tests on a multi-rotationally symmetric body model, which consists of five conical targets. The rotation axes are (0.295, 0.683, 0.669), (-0.512, 0.307, 0.801), (0.426, -0.781, 0.453), (-0.174, -0.592, 0.785), and (0.731, -0.218, 0.644), with axis center coordinates of (56.23, -42.71, 78.19), (-31.85, 69.42, -27.56), (82.14, 15.69, -53.87), (-67.38, -29.54, 41.92), and (12.76, -81.35, 22.69). The first three cones are metal, and the last two are dielectric materials. Figure 6 As shown, a comparison of the magnitudes of unknowns in RCS calculations at different frequencies is presented between the second-order EPA-BoR algorithm and the first-order algorithm; Figure 7 As shown, a comparison of the measured time consumption of the second-order EPA-BoR algorithm and the first-order algorithm for RCS calculation at different frequencies is presented. Under the same computing resource configuration and parameter settings, the proposed algorithm achieves a time saving rate of up to 59.6% compared to the first-order algorithm.

[0167] All numerical experiments in this embodiment were performed on a high-performance computing node equipped with a multi-accelerator architecture. The specific configuration of the node is as follows: the CPU is a Hygon C86 7185 processor (32 cores @ 2.5 GHz), with a memory capacity of 128 GB; it integrates four Hygon DCU-Z100 computing accelerator cards, with an operating frequency of 1600 MHz, and each card is equipped with 32 GB of video memory.

[0168] The core innovation of this invention lies in its excellent parallel scalability and computational efficiency, making it particularly suitable for heterogeneous accelerated computing platforms. The introduction of higher-order line basis functions effectively reduces matrix size and memory usage while ensuring solution accuracy. Based on the DCU architecture, it completes the adaptation of higher-order matrix filling, cyclic reduction, and boundary partitioning, realizing LU decomposition and inversion calculation under the architecture, which significantly reduces memory requirements and communication overhead.

[0169] The above are merely preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle, characterized in that, Includes the following steps: Step 1: Perform geometric modeling and equivalent surface modeling for rotationally symmetric targets: For the cluster of targets to be analyzed, perform independent geometric modeling for each rotationally symmetric target; Based on Huygens' equivalence principle, a closed equivalent surface is constructed outside each rotationally symmetric target; the generatrices of the rotationally symmetric target and the equivalent surface are discretized. Step 2: Construct the EPA-BoR equation system and derive the higher-order U matrix and transmission matrix: Establish integral equations for the modeled rotationally symmetric target and equivalent surface; construct the corresponding EPA-BoR equation system based on the method of moments; derive and establish the higher-order U matrix and transmission matrix based on the formulas for higher-order basis functions and higher-order impedance matrices; substitute the formulas for the higher-order U matrix and transmission matrix into the EPA-BoR equation system to construct the higher-order EPA-BoR equation system. Step 3: Generate each matrix of the high-order EPA-BoR equation system based on the heterogeneous parallel computing strategy: For the established high-order EPA-BoR equation system, use the MPI process to divide it into blocks. Each MPI process independently calculates the corresponding scattering operator and translation operator according to the incident field information. For the high-order matrices contained in the scattering and translation operators, a DCU architecture is used to perform heterogeneous filling calculations, cyclic reduction of the solution results of each thread, and boundary partitioning. The calculated high-order impedance matrix is ​​then decomposed and inverted using a DCU architecture. Step 4: After each MPI process generates its corresponding scattering and translation operators, the electromagnetic scattering solution is completed through inter-process transfer calculation: the preliminary equivalent surface electromagnetic flow solution is obtained through MPI transfer calculation; the error is reduced by using an iterative solution method to obtain the final local solution for each subdomain; Step 5: Calculate the radar cross section of the rotationally symmetric target based on the current solution.

2. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in claim 1, is characterized in that, Step 1 includes: Step 1.1: Model each rotationally symmetric target in the cluster of targets to be analyzed independently. A rotationally symmetric target is a spatial geometric model formed by rotating any continuous curve without repeating points around a straight line. The continuous curve is the generatrix, and the straight line is the axis of rotation. Step 1.2: Based on Huygens' equivalence principle, construct a closed equivalent surface that is coaxial with the target and also rotationally symmetric outside each rotationally symmetric target. The equivalent surface is spaced at a set distance from the target surface. Step 1.3: Discretize the generatrices of the rotationally symmetric target and the equivalent surface. After introducing higher-order line basis functions, set the discretization length to three-tenths of the wavelength corresponding to the incident wave frequency. Step 1.4: Number the points and line segments on the generatrix, and record the vector coordinates of the normal vector of each line segment and the three-dimensional coordinates of the midpoint of the line segment.

3. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in claim 2, is characterized in that... Step 2 includes: Step 2.1, Establishment of integral equation: When the incident wave irradiates the rotationally symmetric target, an induced electromagnetic current is generated on the target surface. The induced electromagnetic current generates a scattered field in the far region. Based on the relationship between the incident wave and the surface electromagnetic current, an integral equation of the electromagnetic field is established. The electromagnetic current on the target surface is divided into two directions: tangential and normal. Step 2.2, Establishment of the EPA-BoR equation set: Based on the principle of the method of moments, the constructed electromagnetic field integral equation is transformed into a matrix equation form; the matrix equations are integrated to construct the EPA-BoR equation set describing the relationship between the equivalent surface scattered current and the incident current. Step 2.3: Based on the formulas for higher-order basis functions and higher-order impedance matrices, derive and establish higher-order U matrices and transmission matrices: The U matrix is ​​a tridiagonal matrix, which is extended to a pentagonal matrix after applying second-order basis functions; the transmission matrix includes the transmission matrix between the equivalent surface and the internal rotationally symmetric target, and the transmission matrix between different equivalent surfaces. The transmission matrix is ​​composed of KL transfer operators. The KL operator form is extended to a higher order and filled into the original transmission matrix to form a higher-order transmission matrix. Step 2.4: Substitute the higher-order impedance matrix and the derived higher-order U matrix and higher-order transmission matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations.

4. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in claim 3, is characterized in that... Step 2.2 is as follows: calculate When there are rotationally symmetric targets =1,…, , =1,…, , , They are the first Scattered current density and magnetic flux density on the equivalent surface of a rotationally symmetric target. , The score is the first Incident current density and magnetic flux density on the equivalent surface of a rotationally symmetric target , They are the first Scattered current density and magnetic flux density on the equivalent surface of a rotationally symmetric target. It is the first Scattering operator for a rotationally symmetric target It is the first The rotationally symmetric target and the first The translation operator between the rotationally symmetric target equivalent surfaces, and the EPA-BoR equations are as follows: The matrix equation is as follows: in, For the first The current coefficient matrix of a rotationally symmetric target is the U matrix. For the first A rotationally symmetric K-operator describing the source field relationship. For the first An L-operator that describes the source field relationship using rotationally symmetric targets. For the first The impedance matrix of a rotationally symmetric target; It is the unit vector normal to the equivalent surface, pointing to the region where the field point is located.

5. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle according to claim 4, characterized in that, Step 2.3 is as follows: (1) Number the second-order basis functions. If the number of segments of the generatrix is ​​N, the second-order basis functions are numbered as the starting number 1. Then the first-order basis functions and the second-order basis functions are numbered alternately, and the second-order basis functions are numbered as the ending number 2N-1. The expression for the second-order basis functions is as follows: in, , These represent the first and second segments of the piecewise basis function one, respectively. Indicates the second basis function; Represents the first in the parametric coordinate system part The tangent vector at the point; Indicates the coordinates of the parameters; The expression for the elements of matrix U is as follows: in, To verify the index of the basis functions, represents the index of the basis functions in the current expansion. Indicates definition at position And the serial number is The test basis functions, Indicates definition at position And the serial number is The source basis functions; Substituting the second-order basis function expression into the element expression of the U matrix, each element of the second-order U matrix is ​​the inner product of basis functions. When the two basis functions do not intersect and their inner product is zero, the U matrix is ​​expanded into a pentagonal matrix, as follows: in, The basis functions are represented sequentially as follows: 1, 2, 3, ..., 2N-3, 2N-2, 2N-1; (2) Based on the different types of fields and flows, the KL transfer operators in the transfer matrix are divided into two categories: electric field The operators relating to the interaction between current and magnetic current are respectively , ; magnetic field The operators relating to the interaction between current and magnetic current are respectively and ; in, , These represent the field point position vector and the source point position vector, respectively. Expressions and Consistent; Expressions and In comparison, free space conductivity Replace with free space permeability ; Based on the method of moments, the matrix form of the KL transfer operator consists of four parts, which are the result of the interaction between the current in the tangential direction and the current in the circumferential direction of the busbar. The expressions of the sub-matrices are as follows: in, Represents electric field With magnetic current The interaction matrix between them Represents electric field With current The interaction matrix between them Represents the first in the matrix OK Column elements; Indicates definition at position And the serial number is The source basis functions; Indicates definition at position And the serial number is The source basis functions; , These represent the surface element vector at the field point and the surface element vector at the source point, respectively. For free space wavenumber, Permeability in free space Let be the Euclidean linear distance between the field point and the source point. For the time-harmonic field scalar Green's function in three-dimensional free space, Angular frequency, The phase constant, Represents the imaginary unit; Represents the gradient operator; , These represent the surface divergence operators acting at the field point and the source point, respectively. Substituting the expression for the second-order basis function into the transfer matrix, the matrix... OK The elements at each column are summed from the calculation results of each subdivision segment. Based on the properties of higher-order line basis functions, only the test basis function is used. and source basis functions The subdivision segment contains non-zero values; since the basis functions span two adjacent subdivision segments, in actual solution, the first... OK The elements in the column are the result of the interaction of four elements from two adjacent subdivision segments.

6. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in claim 5, is characterized in that... Step 3 includes: Step 3.1: For the established high-order EPA-BoR equation set, use the MPI process to divide it into blocks: each MPI process independently reads the geometric and material parameters of the assigned BoR target and classifies the target as a metallic target or a dielectric target, and establishes different MPI task divisions; each MPI process independently calculates the corresponding scattering operator and translation operator based on the incident field information. Step 3.2: For the higher-order matrices contained in the scattering and translation operators, a DCU architecture is used for heterogeneous filling calculations: No inter-process communication is required during the matrix filling stage; each process independently completes the filling of its assigned submatrix. Matrix filling is achieved by iteratively traversing the partitioned segments corresponding to the basis functions and test functions, as detailed below: For the allocation Each test function and For a single MPI process with basis functions, the matrix elements corresponding to the computation are as follows: in, The serial number is the Fourier model number. , This represents the highest order of the Fourier module; Indicates the first The first Fourier mode OK Column matrix elements; Indicates the integral kernel The applied linear operator; This indicates the node number of the Gaussian integral of the test function. This indicates the node number of the Gaussian integral of the basis functions. Indicates the index of the Gaussian integration node of the integration kernel; This represents the total number of nodes in the Gaussian integral of the test function. This represents the total number of nodes in the Gaussian integral of the basis functions. This represents the total number of Gaussian integration nodes in the integration kernel; Indicates the first Gaussian integration weight coefficients corresponding to the Gaussian integration nodes of the test function. Indicates the first Gaussian integration weight coefficients corresponding to the Gaussian integration nodes of the basis functions. Indicates the first Gaussian integration weight coefficients corresponding to each Gaussian integration node of the integration kernel; First, all basis functions required for matrix filling are generated and transferred to the accelerator's global memory. Then, the accelerator maps the integration calculation task to various threads based on the matrix dimensions and integration node information, where the grid dimensions are determined by... Configure the thread block dimension by Configure; Step 3.3: Perform heterogeneous cyclic reduction and boundary partitioning on the solution results of each thread in the DCU: After each thread completes the calculation of its own integral value, iterative updates of the cyclic reduction are performed using a cyclic reduction scheme; after the reduction is completed, the boundary is re-partitioned based on the properties of higher-order basis functions, as follows: After each thread calculates the integral value, the result is reduced to form the complete matrix elements. A cyclic reduction scheme is adopted, and a one-dimensional index is reassigned to the integral value of each thread. The iterative update of the cyclic reduction is completed according to the following formula: Among them, the one-dimensional index in the specification process The three dimensions of a thread block are as follows: ; This refers to the sequence number of the protocol round. Indicates the first The index in the round specification is The integral value, Indicates the first The index in the round specification is The integral value, Indicates the first The index in the round specification is The integral value; Through iterative reduction by pairwise grouping round after round, the final result is achieved... After round reduction, the index is obtained. The integral value of 1 This is the final summation of the integral values ​​of each thread corresponding to the matrix element; Based on the properties of second-order basis functions, when filling matrix elements at the end of reduction accumulation, to ensure that each DCU thread operates completely independently, each MPI process is responsible for calculating the matrix element values ​​on a segment of the busbar; the basis functions of adjacent processes are second-order. They will not intersect; the first segment of piecewise basis function one. The second segment of segmented basis function one, allocated to the previous process. The matrix is ​​assigned to the next process; after each process has calculated its corresponding block matrix, the matrix boundaries of two adjacent processes are overlapped and added together to obtain the result of the corresponding part of the original matrix. Step 3.4: For the calculated high-order impedance matrix, perform LU decomposition and inversion calculation using a DCU architecture: In the direct solver based on LU decomposition, the three core computationally intensive operations, LU decomposition, triangular matrix inversion, and matrix multiplication, are all executed at the DCU end; the heterogeneous computing interface basic linear algebra subroutine library hipBLAS is used to accelerate LU decomposition and matrix multiplication. For the triangular matrix inversion operation in the second step of parallel decomposition, an acceleration operator is designed to speed up the triangular matrix inversion part on the DCU, as follows: Based on a uniform column-strip partitioning strategy with fixed block granularity, two types of triangular matrices are solved using delayed-correction recursive frameworks from right to left and from left to right, respectively; the triangular matrices to be inverted are uniformly denoted as... , The unit lower triangular matrix corresponding to LU decomposition , The non-singular upper triangular matrix corresponding to LU decomposition According to fixed block particle size Will 1-th order matrix Evenly divided into Each column of blocks, the inverse matrix to be solved. Adopted and Completely consistent column-bar block partitioning rules and initial values ​​set to identity matrix of order The unified delay correction recursive formula is as follows: In the formula Number the columns and blocks. Inverse matrix The Each column of bars, For matrix The Middle Each column block corresponds to The inverse matrix of the main diagonal sub-block. for An identity matrix of order 1. For the first Excluding the main diagonal pieces from each column of blocks. The strictly triangular submatrix after that, This is the set of inverse sub-blocks that have been solved in the recursive direction; when When it is a unit lower triangular matrix, The set of column numbers to the right of the current column block. The recursive process is numbered by column block. from Execute from right to left until step 1; when When it is a non-singular upper triangular matrix, The set of column numbers to the left of the current column block. The recursive process is numbered by column block. From 1 to Execute from left to right; At the thread mapping level, a two-dimensional Cartesian topology is used to construct the thread grid and thread blocks: first, construct... Two-dimensional thread block, The number of threads is allocated to the corresponding side length of the thread block, and then a corresponding two-dimensional grid dimension is constructed based on the number of matrix blocks, so that the two-dimensional coordinates of the threads and the row and column indices of the matrix form a one-to-one mapping relationship; in terms of block partitioning strategy design, a two-level block partitioning scheme combining coarse-grained and fine-grained approaches is adopted: the first level is... The coarse-grained columnar block partitioning decomposes the large-scale global matrix into sub-block tasks adapted to the shared memory capacity on the DCU chip; the second level is... Fine-grained micro-block partitioning, dividing each The sub-blocks are further broken down into The micro-computing units are mapped to a single thread for completion, and all recursive calculations within the micro-block are executed in the thread's private registers, adapting to the execution pipeline characteristics of the DCU.

7. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in claim 6, is characterized in that... Step 4 includes: Step 4.1: Solve the initial equivalent surface electromagnetic flow solution through MPI transfer calculation: Each MPI process solves the initial equivalent surface current on the corresponding equivalent surface ES, which will be used as the initial solution of the iterative solver; the target initialization is completed in each MPI process to preserve the locality of electromagnetic interaction; Step 4.2: Use iterative solution method to reduce error and find the final local solution of each subdomain: In the iterative solution stage, the generalized minimum residual (GMRES) method is adopted to update the local solution of each subdomain through scattering operator and translation operator until the relative residual norm converges to below the preset threshold.

8. The high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in claim 7, is characterized in that, Step 5 includes: Step 5.1: After completing the iterative solution, set the equivalent surface electromagnetic currents obtained by convergence in each subdomain as independent Huygens equivalent radiation sources. Step 5.2: For each equivalent radiation source, use the far-field radiation integral formula to calculate the far-field scattered electric field vector of the equivalent radiation source in the specified observation direction; then, superimpose the far-field scattered electric field vectors generated by all equivalent radiation sources in the same observation direction to obtain the total scattered electric field vector of the entire rotationally symmetric target in that observation direction. Step 5.3: Based on the standard definition of radar cross section, calculate the radar cross section value of the target in the observation direction, defined as the ratio of the power density scattered by the target at the receiver per unit solid angle to the incident power density of the incident wave on the target surface. times.

9. A high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation system based on the method of moments and the equivalence principle, characterized in that, This system is used to implement the high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in any one of claims 1 to 8, the system comprising: Geometric modeling and equivalent surface modeling module for rotationally symmetric targets: For the cluster of targets to be analyzed, each rotationally symmetric target is independently geometrically modeled; based on Huygens' equivalence principle, a closed equivalent surface is constructed outside each rotationally symmetric target; the generatrices of the rotationally symmetric targets and the equivalent surface are discretized. The module for constructing higher-order EPA-BoR equations establishes integral equations for the modeled rotationally symmetric target and equivalent surface; constructs the corresponding EPA-BoR equations based on the method of moments; derives and establishes the higher-order U matrix and transmission matrix based on the formulas for higher-order basis functions and higher-order impedance matrices; and substitutes the formulas for the higher-order U matrix and transmission matrix into the EPA-BoR equations to construct the higher-order EPA-BoR equations. The matrix generation modules for the higher-order EPA-BoR equations are as follows: For the established higher-order EPA-BoR equations, the MPI process is used to divide them into blocks. Each MPI process independently calculates the corresponding scattering operator and translation operator based on the incident field information. For the higher-order matrices contained in the scattering operator and translation operator, the DCU architecture is used to perform heterogeneous filling calculations, cyclic reduction of the solution results of each thread, and boundary partitioning. For the calculated higher-order impedance matrix, the DCU architecture is used to perform LU decomposition and inversion calculations. The electromagnetic scattering solution module generates its own corresponding scattering and translation operators in each MPI process, and then completes the electromagnetic scattering solution through inter-process transfer calculation: the preliminary equivalent surface electromagnetic flow solution is obtained through MPI transfer calculation; the error is reduced by using an iterative solution method, and the final local solution of each subdomain is obtained; The radar cross section calculation module calculates the radar cross section of a rotationally symmetric target based on the current solution results.

10. A mobile terminal, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the high-order rotationally symmetric target heterogeneous parallel electromagnetic simulation method based on the method of moments and the equivalence principle as described in any one of claims 1 to 8.

Images