CN118069969A - GPU-based method and device for fast calculation of Green's function in layered media - Google Patents

Abstract

The application provides a method and a device for quickly calculating a hierarchical media green function based on a Graphic Processing Unit (GPU), and relates to the technical field of computational electromagnetism, wherein the method comprises the following steps: initializing a three-dimensional grid of the GPU and the thread number of each thread block; using a calculation task filling matrix of Somoprofil integration (Sommerfeld Integral, SI) of a plurality of parameter points and a plurality of space points contained in the initialized GPU, promoting the numerical integration of SI into matrix products, uniformly distributing calculation tasks of entries of the matrix into each thread block for parallel execution, and obtaining SI calculation results of the plurality of parameter points and the plurality of space points once, wherein the SI calculation results comprise SI head integration results and tail integration results, and the calculation process in each thread block comprises the following steps: and (3) performing matrix product by using a CUDA matrix operation unit Tensor Core, calculating the segment integral of the head and the tail, and accelerating convergence of the segment integral result by using Euler transformation when the tail integral is calculated. The method and the device for calculating the hierarchical medium greens by adopting the scheme realize rapid calculation of the hierarchical medium greens.

Description

GPU-based method and device for fast calculation of Green's function in layered media

Technical Field

The present application relates to the technical field of computational electromagnetics, and in particular to a method and device for fast calculation of Green's functions of layered media based on a GPU.

Background technique

The integral equation method in planar layered media is one of the most successful models in computational electromagnetics and has been widely used in the analysis of microstrip, RF circuits and chips. The difficulty of using the moment method to solve the target electromagnetic response in layered media is mainly the fast calculation of the Green's function. Unlike free space, the Green's function in layered media cannot be written analytically, and the spectral domain Green's function needs to be converted into the spatial domain Green's function through the calculation of the Sommerfeld Integral (SI). However, the high oscillation and slow decay of the Bessel function in the integrand and the singularity of the integral kernel itself make it very difficult to calculate, and its computational efficiency directly affects the filling time of the moment method matrix equation. Therefore, efficiently solving the Sommerfeld integral is the key to speeding up the electromagnetic simulation of planar layered media.

The wide range of applications of Sommerfeld integral continues to attract the attention of researchers. Currently, there are many methods that can accelerate the calculation of Sommerfeld integral. These methods can generally be divided into two categories: closed-form approximation methods and numerical integration methods.

Closed-form approximate method to calculate Sommerfeld integral implementation:

(1) Fitting the spectral domain Green’s function;

(2) By transforming the integral identity into the spatial domain, we can obtain the form of superposition of several spherical waves and/or cylindrical waves;

The most representative closed-form approximation methods are the discrete complex image method (DCIM) and the rational fitting method. Although the closed-form approximation method avoids the calculation of infinite oscillation integrals and greatly reduces the calculation time, the accuracy is uncontrollable, and it is difficult to accurately locate the position of the surface wave pole in layered media.

The implementation scheme for solving the Sommerfeld integral by numerical integration method is as follows:

(1) Determine the integration path;

(2) Direct numerical integration along the integral path;

Numerical integration methods mainly include the Steepest Descent Path (SDP) method, the fast Hankel transform method, and a series of processing methods for the slow decay of the integral kernel. Taking the Steepest Descent Path method as an example, this method mainly deals with the exponential term of the integral kernel, selects the integral path according to the saddle point, and makes the exponential function drop rapidly from the saddle point. The disadvantage of this method is that the integral kernel of the Sommerfeld integral may contain multiple different exponential terms, and the integral path must include each saddle point. When the number of layers is large, the number of saddle points increases rapidly, so it is not suitable for solving general multi-layer structure problems.

Summary of the invention

The present application aims to solve one of the technical problems in the related art at least to some extent.

To this end, the first purpose of this application is to propose a GPU-based fast calculation method for Green's function in layered media, which solves the technical problem that existing methods are difficult to apply in multi-layer structures, realizes the one-time parallel calculation of Sommerfeld integrals of multiple parameters, and can greatly improve the calculation efficiency of Sommerfeld integrals in Green's function in layered media.

The second objective of the present application is to propose a GPU-based fast calculation device for layered medium Green's function.

To achieve the above-mentioned purpose, the first aspect of the present application proposes a GPU-based fast calculation method for layered medium Green's function, including: initializing the three-dimensional grid of the GPU and the number of threads in each thread block; using the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU to fill the matrix, generalize the numerical integral of SI to matrix product, and evenly distribute the calculation tasks of the matrix items to each thread block for parallel execution, so as to obtain SI calculation results of multiple parameter points and multiple spatial points at one time, wherein the SI calculation results include SI head integral results and tail integral results, and the calculation process in each thread block includes: using the CUDA matrix operation unit Tensor Core to perform matrix product, calculate the piecewise integral of the head and tail, and when calculating the tail integral, use Euler transform to accelerate the convergence of the piecewise integral results.

The GPU-based fast calculation method of layered media Green's function in the embodiment of the present application optimizes the calculation architecture by reusing the spectral domain Green's function and Bessel function in the scanning process, converts the Somofi integral into the form of two matrix multiplications, and utilizes the powerful parallel computing capability of the GPU and the special matrix operation unit Tensor Core, as well as the derivation of the Euler transform expression to accelerate the tail integral to realize a parallel solution for the Somofi integral when calculating multiple frequencies or multiple plane layered media parameters at one time.

Optionally, in one embodiment of the present application, the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU are used to fill the matrix, and the numerical integration of SI is generalized to matrix product, including:

Setting the initialization GPU includes parameter points, /> The calculation task of the Somofi integral of the spatial points determines that each thread block calculates the SI of M parameter points and N spatial points;

The SI calculated by each thread block is arranged in an M×N matrix, where each column represents the SI with a different parameter point and each row represents the SI with a different spatial point, so that the numerical integration of SI is generalized to matrix product, and a first matrix and a second matrix are obtained.

Optionally, in one embodiment of the present application, the matrix product is ,/> , the first matrix/> For/> Matrix, the first matrix is composed of the spectral domain Green's function of M parameter points and K integral sampling points, the second matrix/> The term consists of the calculation result of the product of the Bessel function and the integral weight coefficient.

Optionally, in one embodiment of the present application, the Somofi integral includes a head integral and a tail integral, both of which are piecewise integrals, and the head integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the lateral distance between the field point and the source point, /> is the first kind Bessel function, /> is the order of the Bessel function, A is the major axis, /> and/> are weights and samples respectively,/> represents the SI integral result of the i-th sampling point along the elliptical path, and N represents the number of integral sampling points;

The tail integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the first kind Bessel function, /> is the order of the Bessel function, /> is the lateral distance between the field point and the source point, A is the major axis, /> and/> Represent weights and sampling points respectively, L represents the number of sampling points of the sub-intervals of the tail integral interval, N represents the number of sampling points of the sub-integral intervals, /> Represents the calculation result after Euler transformation of the tail integral subinterval.

Optionally, in one embodiment of the present application, the SI tail integral after the Kth recursion during calculation is expressed as:

Where N represents the number of sampling points of the tail integral division subinterval, The coefficients representing the corresponding piecewise integral values, /> Indicates the integral value of each segment of the tail integral.

To achieve the above-mentioned purpose, the second embodiment of the present application proposes a GPU-based layered medium Green's function fast calculation device, including a CPU and a GPU, wherein the CPU includes a memory, wherein:

The CPU is used to initialize the GPU, fill the matrix with the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU, generalize the numerical integral of SI to matrix product, and store the generalized data in the memory;

The GPU is used to evenly distribute the calculation tasks of the matrix items to each thread block for parallel execution, obtain the SI calculation results of multiple parameter points and multiple spatial points at one time, and transmit the integral calculation results to the CPU memory through the PCIe bus. The SI calculation results include the SI head integral result and the tail integral result. The calculation process in each thread block includes:

The CUDA matrix operation unit Tensor Core is used to perform matrix multiplication, calculate the head and tail piecewise integrals, and when calculating the tail integral, the Euler transform is used to accelerate the convergence of the piecewise integral results.

Optionally, in one embodiment of the present application, the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU are used to fill the matrix, and the numerical integration of SI is generalized to matrix product, including:

Setting the initialization GPU includes parameter points, /> The calculation task of the Somofi integral of the spatial points determines that each thread block calculates the SI of M parameter points and N spatial points;

The SI calculated by each thread block is arranged in an M×N matrix, where each column represents the SI with a different parameter point and each row represents the SI with a different spatial point, so that the numerical integration of SI is generalized to matrix product, and a first matrix and a second matrix are obtained.

Optionally, in one embodiment of the present application, the matrix product is ,/> , the first matrix/> For/> Matrix, the first matrix is composed of the spectral domain Green's function of M parameter points and K integral sampling points, the second matrix/> The term consists of the calculation result of the product of the Bessel function and the integral weight coefficient.

Optionally, in one embodiment of the present application, the Somofi integral includes a head integral and a tail integral, both of which are piecewise integrals, and the head integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the lateral distance between the field point and the source point, /> is the first kind Bessel function, /> is the order of the Bessel function, A is the major axis, /> and/> are weights and samples respectively,/> represents the SI integral result of the i-th sampling point along the elliptical path, and N represents the number of integral sampling points;

The tail integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the first kind Bessel function, /> is the order of the Bessel function, /> is the lateral distance between the field point and the source point, A is the major axis, /> and/> Represent weights and sampling points respectively, L represents the number of sampling points of the sub-intervals of the tail integral interval, N represents the number of sampling points of the sub-integral intervals, /> Represents the calculation result after Euler transformation of the tail integral subinterval.

Optionally, in one embodiment of the present application, the implementation method of simplifying the Euler extrapolation method is derived by formula, and the SI tail integral after the Kth recursion during calculation is expressed as:

Where N represents the number of sampling points of the tail integral division subinterval, The coefficients representing the corresponding piecewise integral values, /> Indicates the integral value of each segment of the tail integral.

Additional aspects and advantages of the present application will be given in part in the description below, and in part will become apparent from the description below, or will be learned through the practice of the present application.

BRIEF DESCRIPTION OF THE DRAWINGS

The above and/or additional aspects and advantages of the present application will become apparent and easily understood from the following description of the embodiments in conjunction with the accompanying drawings, in which:

FIG1 is a schematic diagram of a process flow of a GPU-based fast calculation method for layered medium Green's function provided in Embodiment 1 of the present application;

FIG2 is a geometric representation of a planar layered medium according to an embodiment of the present application;

FIG3 is a Sommerfeld numerical integration path diagram of an embodiment of the present application;

FIG4 is an example diagram of using tensor cores to calculate matrix multiplication according to an embodiment of the present application;

FIG5 is a schematic diagram of a one-time parallel calculation scheme of a multi-parameter Sommerfeld integral according to an embodiment of the present application;

FIG6 is a schematic diagram of a three-layered medium model according to an embodiment of the present application;

FIG. 7 is a diagram of the layered Green's function of the three-layer medium at different frequencies in an embodiment of the present application. With/> A first example graph of the magnitudes and relative errors of the components;

FIG8 is a layered Green's function of a three-layer medium at different frequencies in an embodiment of the present application. With/> A second example graph of the magnitude and relative error of the components;

FIG. 9 is a diagram of the layered Green's function of the three-layer medium at different frequencies in an embodiment of the present application. With/> A third example graph of the magnitudes and relative errors of the components;

FIG. 10 is a diagram of the layered Green's function of the three-layer medium at different frequencies in an embodiment of the present application. With/> A fourth example graph of the magnitudes and relative errors of the components;

FIG. 11 is a layered Green's function of a three-layer medium under different dielectric constants in an embodiment of the present application. and/> A first example graph of the magnitudes and relative errors of the components;

FIG. 12 is a layered Green's function of a three-layer medium under different dielectric constants in an embodiment of the present application. and/> A second example graph of the magnitude and relative error of the components;

FIG. 13 is a layered Green's function of a three-layer medium under different dielectric constants according to an embodiment of the present application. and/> A third example graph of the magnitudes and relative errors of the components;

FIG. 14 is a layered Green's function of a three-layer medium under different dielectric constants according to an embodiment of the present application. and/> A fourth example graph of the magnitudes and relative errors of the components;

FIG. 15 is a layered Green's function of a three-layer medium at different medium thicknesses according to an embodiment of the present application. and/> A first example graph of the magnitudes and relative errors of the components;

FIG. 16 is a layered Green's function of a three-layer medium at different medium thicknesses according to an embodiment of the present application. and/> A second example graph of the magnitude and relative error of the components;

FIG. 17 is a layered Green's function of a three-layer medium at different medium thicknesses according to an embodiment of the present application. and/> A third example graph of the magnitudes and relative errors of the components;

FIG. 18 is a layered Green's function of a three-layer medium at different medium thicknesses according to an embodiment of the present application. and/> A fourth example graph of the magnitudes and relative errors of the components;

FIG19 is a schematic diagram of the structure of a GPU-based layered medium Green's function fast calculation device provided in an embodiment of the present application.

Detailed ways

The embodiments of the present application are described in detail below, and examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present application, and should not be construed as limiting the present application.

In related research, many different integral equation forms and Green's functions are used to analyze some specific layered media problems. For example, the electric field integral equation (EFIE), the mixed potential integral equation (MPIE), etc. MPIE has been successfully used in the calculation of planar microstrip antennas and has achieved high calculation accuracy. The most important step in solving MPIE in layered media is the calculation of Green's function. Starting from Maxwell's equations, the various field quantities in the spatial domain are converted to the spectral domain through Fourier transform. Through derivation, it is found that the result has the same form as the transmission line equation. Therefore, the layered medium is equivalent to a transmission line structure. The general expression of the spectral domain Green's function is derived using the transmission line theory. Then, the spatial form of the Green's function in the layered medium can be obtained by calculating the Sommerfeld integral. When using the moment method based on MPIE to analyze the electromagnetic characteristics of layered media, the Green's function is first combined with the mixed potential integral equation, and then the spatial domain Green's function is used. The RWG basis function is selected to expand the current distribution and the Galerkin method is used to match the matrix equation. Finally, the matrix equation is obtained to obtain the current distribution, thereby obtaining the electromagnetic response of the target in the layered medium. The most commonly used MPIE can be expressed as:

in, 、/> denote the magnetic vector potential and electric scalar potential Green's functions, respectively, is a known current source, /> and/> denote the angular frequency and the permittivity and permeability in free space, respectively,/> The matrix form is expressed as:

FIG1 is a flow chart of a method for fast calculation of Green's function of layered media based on a GPU provided in the first embodiment of the present application.

As shown in FIG1 , the GPU-based layered medium Green's function fast calculation method includes the following steps:

Step 101, initializing the three-dimensional grid of the GPU and the number of threads in each thread block;

Step 102, fill the matrix with SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU, generalize the numerical integration of SI to matrix product, and evenly distribute the calculation tasks of the items of the matrix to each thread block for parallel execution, so as to obtain SI calculation results of multiple parameter points and multiple spatial points at one time, wherein the SI calculation results include SI head integration results and tail integration results, and the calculation process in each thread block includes:

The CUDA matrix operation unit Tensor Core is used to perform matrix multiplication, calculate the head and tail piecewise integrals, and when calculating the tail integral, the Euler transform is used to accelerate the convergence of the piecewise integral results.

In some embodiments, multiple value points of other layered media structure parameters contained in the GPU can also be obtained. By filling the matrix with SI calculation tasks of multiple value points and multiple spatial points, this embodiment can achieve rapid scanning of layered media structure parameters.

In some embodiments, the matrix product is calculated using Tensor Core in CUDA, which greatly improves the computational efficiency compared to numerical integration calculations.

In some embodiments, when each thread block performs calculations, the head integral is calculated first and then the tail integral. The Sommerfeld integral result is the sum of the two. The two are serial processes. The calculation processes of the head and tail segmental integrals are similar. The tail integral has one more step than the head integral to calculate the Euler transform of the segmental integral result.

In some embodiments, in the tail integral calculation, an extrapolation algorithm is generally required to accelerate the convergence of the piecewise integral results. Commonly used extrapolation algorithms include Euler transform, average weighted transform, Levin transform, Shanks transform, and the like.

The GPU-based fast calculation method of layered media Green's function in the embodiment of the present application optimizes the calculation architecture by reusing the spectral domain Green's function and Bessel function in the scanning process, converts the Somofi integral into the form of two matrix multiplications, and utilizes the powerful parallel computing capability of the GPU and the special matrix operation unit Tensor Core, as well as the derivation of the Euler transform expression to accelerate the tail integral to realize a parallel solution for the Somofi integral when calculating multiple frequencies or multiple plane layered media parameters at one time.

Optionally, in one embodiment of the present application, the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU are used to fill the matrix, and the numerical integration of SI is generalized to matrix product, including:

Setting the initialization GPU includes parameter points, /> The calculation task of the Somofi integral of the spatial points determines that each thread block calculates the SI of M parameter points and N spatial points;

The SI calculated by each thread block is arranged in an M×N matrix, where each column represents the SI with a different parameter point and each row represents the SI with a different spatial point, so that the numerical integration of SI is generalized to matrix product, and a first matrix and a second matrix are obtained.

Optionally, in one embodiment of the present application, the matrix product is ,/> , the first matrix/> For/> Matrix, the first matrix is composed of the spectral domain Green's function of M parameter points and K integral sampling points, the second matrix/> The term consists of the calculation result of the product of the Bessel function and the integral weight coefficient.

Optionally, in one embodiment of the present application, the Somofi integral includes a head integral and a tail integral, both of which are piecewise integrals, and the head integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the lateral distance between the field point and the source point, /> is the first kind Bessel function, /> is the order of the Bessel function, A is the major axis, /> and/> are weights and samples respectively,/> represents the SI integral result of the i-th sampling point along the elliptical path, and N represents the number of integral sampling points;

The tail integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the first kind Bessel function, /> is the order of the Bessel function, /> is the lateral distance between the field point and the source point, A is the major axis, /> and/> Represent weights and sampling points respectively, L represents the number of sampling points of the sub-intervals of the tail integral interval, N represents the number of sampling points of the sub-integral intervals, /> Represents the calculation result after Euler transformation of the tail integral subinterval.

Optionally, in one embodiment of the present application, the implementation method of the Euler extrapolation method is simplified by formula derivation, and the SI tail integral after the Kth recursion during calculation is expressed as:

Where N represents the number of sampling points of the tail integral division subinterval, The coefficients representing the corresponding piecewise integral values, /> Indicates the integral value of each segment of the tail integral.

The GPU-based layered medium Green's function fast calculation method of the present application is described in detail below through specific embodiments.

Planar layered medium space refers to the discontinuity of the medium only appearing in one direction of the three-dimensional space, and the medium does not change in the other two directions orthogonal to it. It can usually be represented by a layered medium as shown in Figure 2. The medium interface is located at , No./> The relative permittivity and magnetic permeability of the layers are respectively/> The dielectric constant and permeability of the top layer are / > of air or vacuum, assuming the time harmonic factor is/> .

In the spectral domain, the transmission line equation is used to solve the spectral domain Green's function, and then the spatial domain Green's function can be obtained through a two-dimensional inverse Fourier transform. The form of the Sommerfeld integral is expressed as follows:

(1)

in, represents the spectral domain Green's function, /> ,/> are the vertical coordinates of the field point and the source point respectively;/> Based on the source location (/> and/> ) The transverse wave number obtained by transmission line theory; /> is the first kind Bessel function, /> is the order of the Bessel function; is the lateral distance between the field point and the source point;

As shown in Figure 3, the integration path is divided into head integration and Tail integral/> . The head integral path adopts a semi-elliptical path. The major axis/> Generally greater than the maximum value of the wave number in the layered medium, select , where/> . Select /> The values are as follows:

(2)

in, represents the wave number in free space, and the real axis integration path is used to calculate the tail integral. Then the integral expression in equation (1) can be transformed into:

(3)

SI head integral uses the elliptic integral path in the complex plane. In order to ensure the accuracy of the result, the strategy of piecewise integration is adopted. Using the GAUSS-KRONROD integration rule, the head integral can be expressed as:

in, and/> are weight and sampling respectively;/> Indicates the first/> The Sommerfeld integral result of the sampling points along the elliptical path;/> ;

The Sommerfeld tail integral is a semi-infinite integral on the real axis. In order to speed up the convergence and improve the computational efficiency, the classic Euler transformation technique is used to simplify the truncated infinite integral into an L-segment piecewise integral. The piecewise integral can be written as:

(5)

in, and/> Represent weights and sampling points respectively. Segment interval/> is a Bessel function of the first kind/> The interval between two adjacent zero points.

(1) Optimize computing architecture to accelerate

For the Sommerfeld integral shown in equation (1), it can be seen that the integrand is the product of the spectral domain Green's function and the Bessel function. Among them, the spectral domain Green's function is a function related to the frequency and the layered medium parameters, and is independent of the lateral distance of the field source, while the Bessel function is only related to the lateral distance of the field source. This means that in the multi-parameter calculation of the Sommerfeld integral, for the same field source, different frequencies or different layered medium parameters, only the Bessel function needs to be calculated once. Similarly, for the same frequency or the same layered medium parameters, different field sources, only the spectral domain Green's function needs to be calculated once. By reusing the spectral domain Green's function and the Bessel function, the calculation efficiency of the multi-parameter Sommerfeld integral can be greatly improved.

To quantify this speedup, we take a frequency sweep as an example and make the following assumptions:

1) Indicates the calculation of the spectral domain Green's function of a frequency point in a spatial point/> Calculation time.

2) Indicates the calculation of the first kind of Bessel function of a space point/> Calculation time.

3) Indicates the weighted sum of the above calculation results/> time,/> Indicates/> The total number of integral points. Obviously, /> Usually much lower than/> or/> .

calculate parameter points, /> The total time of the Somofi integration of spatial points can be written as/> , using the conventional single parameter point multi-space point SI parallel computing scheme, the total computing time/> can be written as:

(6)

If the intermediate data is used effectively and the architecture optimization method proposed in this application is used, the calculation time becomes

(7)

Therefore, compared with the scheme of single parameter point SI loop calculation, the theoretical speedup ratio is:

(8)

Obviously, the speedup increases with the number of parameter points. and the number of spatial points/> In practical applications, and/> Usually it is a very large number, especially the number of spatial points/> This means that a very high speedup ratio can be achieved by multiplexing the spectral domain Green's function and the Bessel function.

(2) Matrix operations to accelerate integration

In GPU, parameter points, /> The calculation task of the Somofi integral of 𝑀 spatial points is evenly distributed to each thread block. One thread block calculates the Somofi integral of 𝑀 parameter points and 𝑁 spatial points, a total of/> Thread blocks perform integral calculations simultaneously. In order to calculate the integral with/> The spatial distance and /> The SI frequency scan of frequency points can arrange the SI in an M×N matrix, where each column represents the SI of a different parameter point and each row represents the SI of a point with a different spatial distance. By doing so, the numerical integration of SI (including SI head integral and tail integral) can be generalized to matrix product, and the SI calculation results of multiple parameter points and multiple spatial points can be obtained at one time:

(9)

In formula (9), Matrix/> It consists of the spectral domain Green's function of M parameter points and K integral sampling points, and the matrix/> The terms in are composed of the calculation results of the product of Bessel functions and integral weight coefficients.

In modern Nvidia GPUs, there are two hardware units that can be used to perform matrix multiplications (10): CUDA cores and tensor cores. CUDA cores are the basic processing units on the GPU that can perform simple floating-point operations and are optimized for parallel computing workloads. Tensor cores are newer processing units in CUDA that are specifically designed to accelerate tensor operations widely used in deep learning and artificial intelligence applications. Tensor cores perform matrix multiplications or additions more efficiently than CUDA cores for general-purpose parallel workloads. As of CUDA 12.2, Tensor Cores can only perform one at a time. With another /> Therefore, for each thread block of size/> /> With size /> /> The product needs to be multiplied by blocks, as shown in Figure 4. and/> Divide into small matrices, then use tensor cores to calculate the product of these small matrices separately, and then combine them together to get the matrix/> .

(3) Euler extrapolation acceleration

In the calculation of the tail integral of SI, the Euler extrapolation algorithm is needed to accelerate the convergence of the infinite integral. As described in Algorithm 1, the input of the Euler transform is an L-segment integral value sequence stored in the shared memory of the thread block. Specifically,

The input of Algorithm 1 is: Segmented integral value: , the output of Algorithm 1 is: Tail integral value: /> ,

The process of Algorithm 1 is: ; k = 0; while/> , , End; /> ;return/> .

However, for the calculation of multi-parameter Sommerfeld tail integrals, shared memory is difficult to meet the storage space required for Euler transform, and frequent read and write operations of shared memory may cause bank conflicts. In order to solve the above problems, the formula of Euler transform is derived and the numerical calculation is simplified. According to Algorithm 1, the SI tail after the kth recursion can be explicitly written as:

(10)

It conforms to the arrangement rules of Pascal's triangle.

(11)

From formula (10) and formula (11), we can get:

(12)

In this way, the adaptive loop in Algorithm 1 is avoided and the piecewise integration is relieved. The occupancy of shared memory in thread blocks avoids the problem of decreased computing efficiency due to large memory usage.

(4) The one-time parallel calculation scheme of the multi-parameter Sommerfeld integral is shown in Figure 5. The details are as follows:

1) Initialize the thread hierarchy configuration; initialize the GPU's 3D grid to (32,32,1); the number of threads in each thread block is (32,1,1)

2) Fill the matrix and/> ; Matrix/> and/> The calculation of the middle term is evenly divided into a total of 1024 thread blocks, and the calculation tasks in each block are assigned to 32 threads for parallel execution. The matrix is filled and stored in the register,

3) Calculate the head and tail segment integrals by matrix multiplication; use CUDA cores and tensor cores to perform matrix multiplication to obtain SI calculation results for multiple parameter points and multiple spatial points at one time.

4) In the tail integral calculation, the Euler transform is continued to be calculated for the piecewise integral result to accelerate convergence; the head integral calculation directly converts the obtained matrix Stored in global memory;

5) GPU transfers the data integration calculation results to CPU through PCIe bus;

This application is dedicated to high-performance calculation of layered Green's functions in layered media, and proposes a method for fast scanning of layered Green's function parameters. By reusing the integrand in multiple parameters and relying on the efficient parallel computing capability of GPU, accurate and efficient calculation of layered Green's functions for multiple parameter points or multiple layered media parameters is achieved.

This application is to calculate the layered Green's function of the three-layer microstrip structure and/> Taking the component as an example, multi-parameter calculations are performed on the frequency, relative dielectric constant of the layered medium, and the height of the layered medium to demonstrate the beneficial effects. The layered medium parameters of the example model are shown in Figure 6. The computing platform used in this application is NVIDIA RTX 6000 Ada for GPU1, NVIDIA GeForce RTX4090 for GPU2, and Intel Xeon (R) Platinum 8280 for CPU.

(1) Frequency sweep

In the frequency band from 2 GHz to 8 GHz, 64 parameter points are sampled at equal intervals of 0.1 GHz ( =64). At the interface between free space and medium (/> =0,/> =0), lateral distance of source/> 500,000 spatial points are sampled at equal intervals within the range/> =500,000. Using the task division scheme of 𝑀=64,𝑁=96, the calculation results and relative errors are shown in Figures 7, 8, 9, and 10. Figures 7, 8, 9, and 10 only show the calculation results of three parameter points at 2GHz, 4GHz, and 6GHz among the 64 frequencies. It can be seen from Figures 7, 8, 9, and 10 that,/> With/> The relative error is kept at The following demonstrates that the method proposed in this application has good calculation accuracy.

In order to demonstrate the time advantage of the method proposed in this application, the calculation time of the two methods for calculating the above 64 parameter points and 500,000 spatial points is statistically calculated, as shown in Table 1. It can be seen that under the task division scheme of M=64 and N=96, the method proposed in this application achieves a 1914-fold acceleration in the head integration and a 1226-fold acceleration in the tail integration compared with the conventional method using OpenMP parallelization, shortening the total time from 24960.62 seconds to less than 14.84 seconds, achieving an acceleration of more than 1600 times.

Table 1 Calculation of three-layer media using different methods Time (M=64,N=96)

(2) Relative dielectric constant scanning

At a frequency of 8 GHz and a medium thickness of 0.254 mm, 64 relative dielectric constants were sampled at intervals of 0.1 in the range of 3.2-9.6. =0,/> =0), the lateral distance of the field source 500,000 spatial points are sampled at equal intervals within the range (/> =500,000), calculate the corresponding /> as well as The components, results and relative errors are shown in Figures 11, 12, 13 and 14. Figures 11, 12, 13 and 14 only show the calculation results when the relative dielectric constant is 3.6, 6.6 and 9.6. It can be found that the relative errors of the method proposed in this application are all within It is proved that the method proposed in this application is still applicable to the scanning of dielectric constant, and the calculation time here is basically the same as that of the multi-parameter point method.

(3) Layered media height scanning

At a frequency of 8 GHz and a relative dielectric constant of 9.6, 64 medium thicknesses were sampled at 0.02 mm intervals within the range of 0.254 mm to 1.534 mm. =0,/> =0), the lateral distance of the field source 500,000 spatial points are sampled at equal intervals within the range (/> =500,000), calculate the corresponding /> as well as The results and relative errors are shown in Figures 15, 16, 17 and 18. Figures 15, 16, 17 and 18 only show the calculation results when the medium thickness is 0.254mm, 0.854mm and 1.454mm. It can be seen that the relative error is still maintained at It is proved that the method proposed in this application is still applicable to the simultaneous calculation of the thickness of multiple layered media, and the calculation time here is basically consistent with the multi-parameter point method.

Numerical experiments show that the relative error of the method proposed in this application is kept within Compared with the single parameter point method accelerated by OpenMP in high-end CPU, the time in the three-layer medium structure calculation example was shortened from 24960.62 seconds to less than 14.84 seconds, achieving a 1682-fold acceleration. The method proposed in this application will have a significant acceleration effect on the simulation and parameter optimization of microstrip circuits and microwave integrated circuits.

In order to implement the above embodiments, the present application also proposes a GPU-based layered medium Green's function fast calculation device.

FIG19 is a schematic diagram of the structure of a GPU-based layered medium Green's function fast calculation device provided in an embodiment of the present application.

As shown in FIG19 , the GPU-based layered medium Green's function fast calculation device includes a CPU and a GPU, wherein the CPU includes a memory, wherein:

The CPU is used to initialize the GPU, fill the matrix with the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU, generalize the numerical integral of SI to matrix product, and store the generalized data in the memory;

The GPU is used to evenly distribute the calculation tasks of the matrix items to each thread block for parallel execution, obtain the SI calculation results of multiple parameter points and multiple spatial points at one time, and transmit the integral calculation results to the CPU memory through the PCIe bus. The SI calculation results include the SI head integral result and the tail integral result. The calculation process in each thread block includes:

The CUDA matrix operation unit Tensor Core is used to perform matrix multiplication, calculate the head and tail piecewise integrals, and when calculating the tail integral, the Euler transform is used to accelerate the convergence of the piecewise integral results.

Optionally, in one embodiment of the present application, the SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU are used to fill the matrix, and the numerical integration of SI is generalized to matrix product, including:

Setting the initialization GPU includes parameter points, /> The calculation task of the Somofi integral of the spatial points determines that each thread block calculates the SI of M parameter points and N spatial points;

The SI calculated by each thread block is arranged in an M×N matrix, where each column represents the SI with a different parameter point and each row represents the SI with a different spatial point, so that the numerical integration of SI is generalized to matrix product, and a first matrix and a second matrix are obtained.

Optionally, in one embodiment of the present application, the matrix product is ,/> , the first matrix/> For/> Matrix, the first matrix is composed of the spectral domain Green's function of M parameter points and K integral sampling points, the second matrix/> The term consists of the calculation result of the product of the Bessel function and the integral weight coefficient.

Optionally, in one embodiment of the present application, the Somofi integral includes a head integral and a tail integral, both of which are piecewise integrals, and the head integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the lateral distance between the field point and the source point, /> is the first kind Bessel function, /> is the order of the Bessel function, A is the major axis, /> and/> are weights and samples respectively,/> represents the SI integral result of the i-th sampling point along the elliptical path, and N represents the number of integral sampling points;

The tail integral is expressed as:

in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the first kind Bessel function, /> is the order of the Bessel function, /> is the lateral distance between the field point and the source point, A is the major axis, /> and/> Represent weights and sampling points respectively, L represents the number of sampling points of the sub-intervals of the tail integral interval, N represents the number of sampling points of the sub-integral intervals, /> Represents the calculation result after Euler transformation of the tail integral subinterval.

Optionally, in one embodiment of the present application, the implementation method of simplifying the Euler extrapolation method is derived by formula, and the SI tail integral after the Kth recursion during calculation is expressed as:

Where N represents the number of sampling points of the tail integral division subinterval, The coefficients representing the corresponding piecewise integral values, /> Indicates the integral value of each segment of the tail integral.

It should be noted that the above explanation of the embodiment of the GPU-based layered medium Green's function fast calculation method is also applicable to the GPU-based layered medium Green's function fast calculation device of this embodiment, and will not be repeated here.

In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example" or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, unless they are contradictory.

In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In the description of this application, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

Any process or method description in a flowchart or otherwise described herein may be understood to represent a module, fragment or portion of code comprising one or more executable instructions for implementing the steps of a custom logical function or process, and the scope of the preferred embodiments of the present application includes alternative implementations in which functions may not be performed in the order shown or discussed, including performing functions in a substantially simultaneous manner or in reverse order depending on the functions involved, which should be understood by technicians in the technical field to which the embodiments of the present application belong.

The logic and/or steps represented in the flowchart or otherwise described herein, for example, can be considered as an ordered list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by an instruction execution system, device or apparatus (such as a computer-based system, a system including a processor, or other system that can fetch instructions from an instruction execution system, device or apparatus and execute instructions), or in combination with these instruction execution systems, devices or apparatuses. For the purposes of this specification, "computer-readable medium" can be any device that can contain, store, communicate, propagate or transmit a program for use by an instruction execution system, device or apparatus, or in combination with these instruction execution systems, devices or apparatuses. More specific examples (non-exhaustive list) of computer-readable media include the following: an electrical connection with one or more wires (electronic device), a portable computer disk box (magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable and programmable read-only memory (EPROM or flash memory), an optical fiber device, and a portable compact disk read-only memory (CDROM). In addition, the computer-readable medium may even be paper or other suitable medium on which the program is printed, since the program may be obtained electronically, for example, by optically scanning the paper or other medium and then editing, interpreting or processing in other suitable ways if necessary, and then stored in a computer memory.

It should be understood that the various parts of the present application can be implemented by hardware, software, firmware or a combination thereof. In the above-mentioned embodiments, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, it can be implemented by any one of the following technologies known in the art or their combination: a discrete logic circuit having a logic gate circuit for implementing a logic function for a data signal, a dedicated integrated circuit having a suitable combination of logic gate circuits, a programmable gate array (PGA), a field programmable gate array (FPGA), etc.

A person skilled in the art may understand that all or part of the steps in the method for implementing the above-mentioned embodiment may be completed by instructing related hardware through a program, and the program may be stored in a computer-readable storage medium, which, when executed, includes one or a combination of the steps of the method embodiment.

In addition, each functional unit in each embodiment of the present application may be integrated into a processing module, or each unit may exist physically separately, or two or more units may be integrated into one module. The above-mentioned integrated module may be implemented in the form of hardware or in the form of a software functional module. If the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it may also be stored in a computer-readable storage medium.

The storage medium mentioned above may be a read-only memory, a disk or an optical disk, etc. Although the embodiments of the present application have been shown and described above, it can be understood that the above embodiments are exemplary and cannot be understood as limiting the present application. A person of ordinary skill in the art may change, modify, replace and modify the above embodiments within the scope of the present application.

Claims (10)

1. A GPU-based fast calculation method for layered medium Green's function, characterized in that the method is applied to electromagnetic simulation of radio frequency integrated circuits or devices in planar layered media, and the method comprises the following steps: Initialize the GPU's 3D grid and the number of threads per thread block; The SI calculation tasks of multiple parameter points and multiple spatial points contained in the initialized GPU are used to fill the matrix, the numerical integration of SI is generalized to matrix product, and the calculation tasks of the items of the matrix are evenly distributed to each thread block for parallel execution, so as to obtain the SI calculation results of multiple parameter points and multiple spatial points at one time, wherein the SI calculation results include the SI head integration result and the tail integration result, and the calculation process in each thread block includes: The CUDA matrix operation unit Tensor Core is used to perform matrix multiplication, calculate the head and tail piecewise integrals, and when calculating the tail integral, the Euler transform is used to accelerate the convergence of the piecewise integral results. 2. The GPU-based layered medium Green's function fast calculation method according to claim 1, characterized in that the SI calculation task filling matrix using the multiple parameter points and multiple spatial points contained in the initialized GPU, and generalizing the numerical integration of SI to matrix product, comprises: Setting the initialization GPU includes parameter points, /> The calculation task of the Somofi integral of the spatial points determines that each thread block calculates the SI of M parameter points and N spatial points; The SI calculated by each thread block is arranged in an M×N matrix, where each column represents the SI with a different parameter point and each row represents the SI with a different spatial point, so that the numerical integration of SI is generalized to matrix product, and a first matrix and a second matrix are obtained. 3. The GPU-based fast calculation method for layered medium Green's function according to claim 2, wherein the matrix product is ,/> , the first matrix/> For/> Matrix, the first matrix is composed of spectral domain Green's functions of M parameter points and K integral sampling points, the second matrix/> The term consists of the calculation result of the product of the Bessel function and the integral weight coefficient. 4. The GPU-based layered medium Green's function fast calculation method according to claim 1, characterized in that the Somofi integral includes a head integral and a tail integral, both of which are piecewise integrals, and the head integral is expressed as: in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the lateral distance between the field point and the source point, /> is the first kind Bessel function, /> is the order of the Bessel function, A is the major axis, /> and/> are weights and samples respectively,/> represents the SI integral result of the i-th sampling point along the elliptical path, and N represents the number of integral sampling points; The tail integral is expressed as: in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the first kind Bessel function, /> is the order of the Bessel function, /> is the lateral distance between the field point and the source point, A is the major axis, /> and/> Represent weights and sampling points respectively, L represents the number of sampling points of the sub-intervals of the tail integral interval, N represents the number of sampling points of the sub-integral intervals, /> Represents the calculation result after Euler transformation of the tail integral subinterval. 5. The GPU-based fast calculation method for layered medium Green's function as claimed in claim 1, characterized in that the implementation method of simplifying the Euler extrapolation method is derived by formula, and the SI tail integral after the K-th recursion during calculation is expressed as: Where N represents the number of sampling points of the tail integral division subinterval, The coefficients representing the corresponding piecewise integral values, /> Indicates the integral value of each segment of the tail integral. 6. A GPU-based device for fast calculation of Green's function in layered media, characterized in that the device is applied to electromagnetic simulation of radio frequency integrated circuits or devices in planar layered media, the device comprises a CPU and a GPU, the CPU comprises a memory, wherein: The CPU is used to initialize the GPU, fill the matrix with the calculation tasks of SI of multiple parameter points and multiple spatial points contained in the initialized GPU, generalize the numerical integral of SI to matrix product, and store the generalized data in the memory; The GPU is used to evenly distribute the calculation tasks of the matrix items to each thread block for parallel execution, obtain SI calculation results of multiple parameter points and multiple spatial points at one time, and transmit the integral calculation results to the memory of the CPU through the PCIe bus, wherein the SI calculation results include SI head integral results and tail integral results, and the calculation process in each thread block includes: The CUDA matrix operation unit Tensor Core is used to perform matrix multiplication, calculate the head and tail piecewise integrals, and when calculating the tail integral, the Euler transform is used to accelerate the convergence of the piecewise integral results. 7. The GPU-based layered medium Green's function fast calculation device according to claim 6, characterized in that the SI calculation task filling matrix using the multiple parameter points and multiple spatial points contained in the initialized GPU and generalizing the numerical integration of SI to matrix product comprises: Setting the initialization GPU includes parameter points, /> The calculation task of the Somofi integral of the spatial points determines that each thread block calculates the SI of M parameter points and N spatial points; The SI calculated by each thread block is arranged in an M×N matrix, where each column represents the SI with a different parameter point and each row represents the SI with a different spatial point, so that the numerical integration of SI is generalized to matrix product, and a first matrix and a second matrix are obtained. 8. The GPU-based layered medium Green's function fast calculation device according to claim 7, wherein the matrix product is ,/> , the first matrix/> For/> Matrix, the first matrix is composed of spectral domain Green's functions of M parameter points and K integral sampling points, the second matrix/> The term consists of the calculation result of the product of the Bessel function and the integral weight coefficient. 9. The GPU-based layered medium Green's function fast calculation device according to claim 6, characterized in that the Somofi integral includes a head integral and a tail integral, both of which are piecewise integrals, and the head integral is expressed as: in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the lateral distance between the field point and the source point, /> is the first kind Bessel function, /> is the order of the Bessel function, A is the major axis, /> and/> are weights and samples respectively,/> represents the SI integral result of the i-th sampling point along the elliptical path, and N represents the number of integral sampling points; The tail integral is expressed as: in, represents the spectral domain Green's function, /> 、/> are the vertical coordinates of the field point and the source point, respectively,/> is the transverse wave number obtained by transmission line theory based on the source position,/> is the first kind Bessel function, /> is the order of the Bessel function, /> is the lateral distance between the field point and the source point, A is the major axis, /> and/> Represent weights and sampling points respectively, L represents the number of sampling points of the sub-intervals of the tail integral interval, N represents the number of sampling points of the sub-integral intervals, /> Represents the calculation result after Euler transformation of the tail integral subinterval. 10. The GPU-based layered medium Green's function fast calculation device according to claim 6, characterized in that the implementation method of simplifying the Euler extrapolation method is derived by formula, and the SI tail integral after the K-th recursion during calculation is expressed as: Where N represents the number of sampling points of the tail integral division subinterval, The coefficients representing the corresponding piecewise integral values, /> Indicates the integral value of each segment of the tail integral.