Adaptive swept-frequency electromagnetic simulation method based on accelerated GMRES solver

Abstract

The present invention discloses an adaptive frequency sweeping electromagnetic simulation method based on an accelerated GMRES solver. In electromagnetic simulation calculations, for large-scale and complex examples, the number of RWG basis functions formed after decomposition is as high as tens of thousands, resulting in a very large matrix equation. The present invention proposes an adaptive frequency sweeping electromagnetic simulation method based on an accelerated GMRES solver, which utilizes the characteristic that the current changes smoothly at adjacent frequency points in electromagnetic calculations, adopts an interpolation method to estimate the current value of the unknown point, and then uses the estimated current value as the initial value of the equation for solving the problem. This can make the initial current value closer to the actual solved current value, thereby significantly reducing the size of the initial residual, effectively improving the problem of large errors in the estimated current value by the solver, and thus reducing the number of iteration steps at the current frequency point, and combining it with the adaptive frequency sweeping method in the moment method simulation calculation, effectively improving the simulation efficiency while ensuring the accuracy of the results.

Description

Technical Field

[0001] The present invention belongs to the technical field of radio frequency circuit engineering, and relates to an adaptive frequency sweeping electromagnetic simulation method based on an accelerated GMRES (generalized minimum residual method) solver. Background Art

[0002] With the continuous development of computational electromagnetics and simulation technology, the "electrically large" size of the problems to be solved in current scientific and engineering tasks is getting larger and larger. Since the number of unknowns formed by network decomposition of complex "electrically large" size is too large, and the size of the equation group under the moment method is affected by the number of RWG basis functions formed after network decomposition, the scale of the formed equation group is very large. At this time, if you continue to choose the direct solution method, it is easy to be limited by computer hardware and resources. In the case of limited memory resources, the equation group data will not be able to be stored, resulting in program crash and inability to solve the equation group. An iterative method should be used for solution. Although the iterative algorithm can reduce the time complexity, the large gap between its initialization current value and the actual current value of the simulation leads to too high an iteration step number, so that the simulation time is too slow. The present invention proposes an adaptive sweep frequency electromagnetic simulation method based on an accelerated GMRES solver. The method uses the idea that the current changes smoothly between adjacent frequency points to pre-interpolate the current value of the current at the current simulation frequency point. The predicted current value is close to the actual value to a certain extent, which effectively improves the problem of large error in the solver's estimated current value, and is combined with the adaptive sweep frequency method in the moment method to improve the simulation efficiency. Summary of the Invention

[0003] The purpose of the present invention is to address the deficiencies of the prior art and provide an adaptive swept-frequency electromagnetic simulation method based on an accelerated GMRES solver, which has fast calculation and high accuracy.

[0004] In a first aspect, the present invention provides an adaptive swept-frequency electromagnetic simulation method based on an accelerated GMRES solver, the method comprising the following steps:

[0005] Step (1) Consider an ideal conductor in a uniform space with a dielectric constant of ε and a magnetic permeability of μ. When it is affected by the incident field in any direction, it will generate a corresponding induced electromagnetic current and an outward scattered field. Assume that the incident electric field of the incident plane electromagnetic wave is E inc , the magnetic field is H inc , the scattered electric field generated by irradiating the surface of an ideal conductor is E sac , the scattered magnetic field is H sac , inside an ideal conductor, the internal electric field E in and the internal magnetic field H in All are zero;

[0006] Step (2), assuming the simulation frequency band is set to [f1,f n], according to the bisection interpolation theory, we know that there are three frequency points selected in the first round, namely f1, f n and f mid , f1 is the starting frequency point of the simulation, f n is the end frequency point of the simulation, f mid is the middle value between the start and end frequency points, f mid =(f1+f n ) / 2;

[0007] Step (3), using the moment method to solve the electromagnetic field inside the conductor at the currently selected frequency point;

[0008] According to the electric field integral equation:

[0009]

[0010] In the formula represents the normal unit vector outside the conductor, η represents the wave impedance of the medium space, L represents the linear operator, J represents the current density, E inc represents the incident electric field;

[0011] The RWG basis function is defined as being able to simulate the electromagnetic flow distribution on the surface of an object of arbitrary shape when there are common edges between two adjacent triangles;

[0012] The definition of the RWG basis function is given as:

[0013]

[0014] where l n is the nth pair of adjacent triangles T n + and The length of the common side, r is any point on the triangle, From triangle T n + Vertex V n + The vector pointing to point r, From the triangle Vertex The vector pointing to point r, and They are adjacent triangles T n + and The area of ​​​​f n (r) represents the nth RWG basis function in the field region;

[0015] By expanding the current in equation (1) using the RWG basis function:

[0016]

[0017] Among them I n is the nth current expansion coefficient to be solved, and N is the number of unknowns formed after the network is partitioned;

[0018] Substituting formula (3) into formula (1) yields the discretized electric field integral equation:

[0019]

[0020] Where L(f n (r)) represents the RWG basis function f n (r) linear operator equation;

[0021] The Galerkin method is used to select the RWG basis function as the test function to test the electric field integral equation (4):

[0022]

[0023] where f n (r) represents the nth RWG basis function in the field region; f m (r') represents the m-th RWG basis function at point r' in the source region;

[0024] Since the RWG basis function is along the surface tangent direction, formula (5) can be rewritten as follows:

[0025]

[0026] Write formula (6) in the form of a matrix equation:

[0027] ZI=V Formula (7)

[0028] Where Z is the N×N moment method impedance matrix, V and I are N×1 column vectors, representing the excitation vector and current vector in the moment method, respectively;

[0029] The elements of the mth row and nth column of the matrix Z and the mth row of the matrix V are expressed as follows:

[0030]

[0031]

[0032]

[0033] Where j represents the imaginary unit, f represents the frequency point; RWG basis function f m The source triangle pair where (r') is located is RWG basis function f n The field triangle where (r) is located is G(R) is the Green's function in uniform space, R = |r-r'| represents the distance from any point r of the field basis function to any point r' of the source basis function; represents the gradient operator, represents the divergence operator, ds represents the surface vector element of the field region, and ds' represents the surface vector element of the source region;

[0034] Step (4), using the accelerated GMRES (generalized minimum coparameter) solver to iteratively solve the matrix equation (7) to obtain the current value at the simulation frequency point;

[0035] 4-1 Calculate the current value at the current frequency point

[0036] 4-1-1 For a matrix equation (7), initialize I i =I0, I0 is initialized to 0, the number of iterations i = 0;

[0037] 4-1-2 According to formula (7), we can know formula (11), and then get the residual r i ;

[0038] r i =V-ZI i Formula (11)

[0039] 4-1-3 Determine the residual r i Is it less than the artificially set threshold k? If so, it is considered that the I i This is the current value you are looking for, otherwise go to step 4-1-4;

[0040] 4-1-4 will be the residual r i The orthogonalized result is used as the basis vector of the initial Krylov subspace, and then the Arnoldi iteration is performed to obtain an approximate solution vector I i+1 ; Then return to step 4-1-2 and update the number of iterations i=i+1;

[0041] 4-2 Calculate the current value at the third frequency point

[0042] 4-2-1 Calculate step (2) f1 and f according to step 4-1 n and f mid The current value of any two frequency points in the calculation is used to calculate the estimated current value of the third frequency point using the interpolation algorithm. i Initialize to the estimated current value at the third frequency point, and initialize the number of iterations i = 0;

[0043] 4-1-2 According to formula (11), the residual r is calculated i ;

[0044] 4-1-3 Determine the residual ri If it is less than 1, then the estimated current value at the current frequency point is considered to be an effective value, otherwise I i Reset to 0 and repeat steps 4-1-2 to 4-1-4;

[0045] Step (5): Simulation accuracy evaluation

[0046] For step (2) f1, f n and f mid The current value results are used for convergence judgment. If the convergence conditions are met, the simulation calculation is stopped. If the convergence conditions are not met, the next round of increasing the sampling frequency is carried out using the bisection interpolation theory.

[0047] rate point, repeat steps (3)-(5);

[0048] The convergence criterion is the error formula Where S(f) represents the theoretical value of frequency response, Represents the frequency response fitting value; if err < threshold, convergence is achieved, otherwise convergence is not achieved;

[0049] Step (6): The electromagnetic simulation parameters can be obtained based on the current value obtained by the above calculation to realize electromagnetic simulation.

[0050] In a second aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed in a computer, causes the computer to execute the method described above.

[0051] In a third aspect, the present invention provides a computing device comprising a memory and a processor, wherein the memory stores executable code, and when the processor executes the executable code, the method described is implemented.

[0052] Based on the characteristic of smooth current change, the present invention uses the current value of the completed frequency point to predict the current value of the frequency point where the simulation has not started, and then uses the predicted value as the initial value. This can effectively reduce the problem of a large gap between the initial current value of the GMRES solver and the actual current value, thereby reducing the number of subsequent iterative steps and avoiding slow simulation time due to incorrect initial current value of the solver. At the same time, the solver is combined with the adaptive frequency sweep method in the moment method simulation to effectively improve the simulation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 This is a flow chart of the a priori initial value method based on accelerated GMRES solution of the present invention;

[0054] Figure 2 This is a structural diagram of the transformer of the experimental example of the present invention;

[0055] Figure 3 This is a layer structure diagram of the transformer of the experimental example of the present invention;

[0056] Figure 4 This is the experimental example S of the present invention 11 Parameter comparison chart;

[0057] Figure 5 This is a comparison chart of the number of iterative steps of the present invention. DETAILED DESCRIPTION

[0058] The present technology will be further described in detail below with reference to the accompanying drawings and implementation cases.

[0059] like Figure 1 As shown in FIG, it is a flow chart of an adaptive frequency sweep electromagnetic simulation method based on an accelerated GMRES solver of the present invention. An example of a transformer is given. Unlike traditional linear windings, the windings of the transformer are distributed in a spiral form. This model can make the structure more compact and achieve a larger inductance value in a limited space. The present invention analyzes from 1 GHz to 20 GHz. The traditional analysis method is to calculate at an average interval of 0.1 GHz. The structure diagram of the transformer is shown in FIG. Figure 2 As shown, its layer structure is as Figure 3 shown.

[0060] Step (1) Consider an ideal conductor in a uniform space with a dielectric constant of ε and a magnetic permeability of μ. When it is affected by the incident field in any direction, it will generate a corresponding induced electromagnetic current and an outward scattered field. Assume that the incident electric field of the incident plane electromagnetic wave is E inc , the magnetic field is H inc , the scattered electric field generated by irradiating the surface of an ideal conductor is E sac , the scattered magnetic field is H sac , inside an ideal conductor, the internal electric field E in and the internal magnetic field H in All are zero;

[0061] Step (2), assuming the simulation frequency band is set to [f1,f n ], according to the bisection interpolation theory, we know that there are three frequency points selected in the first round, namely f1, f n and f mid , f1 is the starting frequency point of the simulation, f n is the end frequency point of the simulation, f mid is the middle value between the start and end frequency points, f mid =(f1+f n ) / 2;

[0062] The first three frequency points are selected as: f1=1GHz, f mid =10.27GHz,f n=20GHz.

[0063] Step (3): Use moment method to solve the internal electromagnetic field of the conductor at the currently selected frequency point

[0064] According to the electric field integral equation:

[0065]

[0066] In the formula represents the normal unit vector outside the conductor, η represents the wave impedance of the medium space, L represents the linear operator, J represents the current density, E inc represents the incident electric field;

[0067] The RWG basis function is defined as being able to simulate the electromagnetic flow distribution on the surface of an object of arbitrary shape when there are common edges between two adjacent triangles;

[0068] The definition of the RWG basis function is given as:

[0069]

[0070] where l n is the nth pair of adjacent triangles T n + and The length of the common side, r is any point on the triangle, From triangle T n + Vertex V n + The vector pointing to point r, From the triangle Vertex V n - The vector pointing to point r, and They are adjacent triangles T n + and The area of ​​​​f n (r) represents the nth RWG basis function in the field region;

[0071] By expanding the current in equation (1) using the RWG basis function:

[0072]

[0073] Among them I n is the nth current expansion coefficient to be solved, and N is the number of unknowns formed after the network is partitioned;

[0074] Substituting formula (3) into formula (1) yields the discretized electric field integral equation:

[0075]

[0076] Where L(f n (r)) represents the RWG basis function f n (r) linear operator equation;

[0077] The Galerkin method is used to select the RWG basis function as the test function to test the electric field integral equation (4):

[0078]

[0079] where f n (r) represents the nth RWG basis function in the field region; f m (r') represents the m-th RWG basis function at point r' in the source region;

[0080] Since the RWG basis function is along the surface tangent direction, formula (5) can be rewritten as follows:

[0081]

[0082] Write formula (6) in the form of a matrix equation:

[0083] ZI=V Formula (7)

[0084] Where Z is the N×N moment method impedance matrix, V and I are N×1 column vectors, representing the excitation vector and current vector in the moment method, respectively;

[0085] The elements of the mth row and nth column of the matrix Z are expressed as follows:

[0086]

[0087]

[0088]

[0089] Where j represents the imaginary unit, f represents the frequency point; RWG basis function f m The source triangle pair where (r') is located is RWG basis function f n The field triangle where (r) is located is G(R) is the Green's function in uniform space, R = |r-r'| represents the distance from any point r of the field basis function to any point r' of the source basis function; represents the gradient operator, represents the divergence operator, ds represents the surface vector element of the field region, and ds' represents the surface vector element of the source region;

[0090] Step (4), using the accelerated GMRES solver to iteratively solve the matrix equation (7) to obtain the current value at the simulation frequency point;

[0091] The accelerated GMRES solver is used to simulate the matrix equation (7) formed by the moment method discretization of the frequency points selected in step (2) to obtain the current value. For any frequency point, the initialization operation must be performed first;

[0092] 4-1 Calculate the current value at the current frequency point

[0093] 4-1-1 For a matrix equation (7), initialize I i =I0, I0 is initialized to 0, the number of iterations i = 0;

[0094] 4-1-2 According to formula (7), we can know formula (11), and then get the residual r i ;

[0095] r i =V-ZI i Formula (11)

[0096] 4-1-3 Determine the residual r i Is it less than the artificially set threshold k? If so, it is considered that the I i This is the current value you are looking for, otherwise go to step 4-1-4;

[0097] 4-1-4 will be the residual r i The orthogonalized result is used as the basis vector of the initial Krylov subspace, and then the Arnoldi iteration is performed to obtain an approximate solution vector I i+1 ; Then return to step 4-1-2 and update the number of iterations i=i+1;

[0098] 4-2 Calculate the current value at the third frequency point

[0099] 4-2-1 Calculate step (2) f1 and f according to step 4-1 n and f mid The current value of any two frequency points in the calculation is used to calculate the estimated current value of the third frequency point using the interpolation algorithm. i Initialize to the estimated current value at the third frequency point, and initialize the number of iterations i = 0;

[0100] For RF microwave circuits, the current distribution varies continuously at different frequencies and across different conductor planes, and the current variation between adjacent frequency points is typically very small. Since the parameter to be calculated in the method of moments is the current value, the current at the unsolved frequency point can be estimated using interpolation based on the current data at known frequency points as the initial current value, and then an iterative solution can be performed.

[0101] 4-1-2 According to formula (11), the residual r is calculated i ;

[0102] 4-1-3 Determine the residual r i If it is less than 1, then the estimated current value at the current frequency point is considered to be an effective value, otherwise I i Reset to 0 and repeat steps 4-1-2 to 4-1-4;

[0103] Step (5): Simulation accuracy evaluation

[0104] For step (2) f1, f n and f mid The current value result is used to make a convergence judgment. If the convergence condition is met, the simulation calculation is stopped. If the convergence condition is not met, the bisection interpolation theory is used to add new sampling frequency points for the next round, and steps (3)-(5) are repeated.

[0105] The convergence criterion is the error formula Where S(f) represents the theoretical value of frequency response, Represents the frequency response fitting value; if err < threshold, convergence is achieved, otherwise convergence is not achieved;

[0106] The new sampling frequency point selection rule is to divide the frequency band into two segments, namely [f1,f mid ] and [f mid ,f n ], in the first half, three sampling points f1, f mid and f mid / 2 , are the starting and ending endpoints and the middle value of the newly divided frequency range respectively. Similarly, the selection rules for the second half are the same as those for the first half;

[0107] Step (6): The electromagnetic simulation S parameters can be obtained based on the current values ​​obtained by the above calculation to realize electromagnetic simulation.

[0108] Step (7), in order to reduce the influence of the Runge phenomenon caused by interpolation, when there are more than a manually set threshold a of frequency points on the left and right sides of the predicted frequency point, the current value of the frequency point farthest from the predicted frequency point is discarded to ensure that the number of current values ​​on both sides is within the threshold a;

[0109] Step (8), loop judgment accuracy evaluation until the entire judgment simulation is completed; step (5) is performed on each set of sampled frequency points until all frequency bands converge;

[0110] like Figure 4 As shown, the S of the accelerated GMRES solver and the original solver 11The parameter result comparison chart shows that the two curves have the same results. Accelerating the GMRES solver does not affect the accuracy in the adaptive frequency sweep method.

[0111] like Figure 5 As shown in Figure 3, the comparison of the number of iteration steps when the accelerated GMRES solver is used, the residual data is calculated every 20 steps in this experiment. Except for the 20 GHz point, the number of steps at other points is reduced, and the efficiency is significantly improved.

[0112] The above embodiments are for demonstration only, and the present invention is not limited to the above embodiments. As long as the requirements of the present invention are met, any similar implementation is included in the protection scope of the present invention.

Claims

1. Adaptive frequency sweep electromagnetic simulation method based on accelerated GMRES solver, characterized by The method comprises the following steps: Step (1) Consider an ideal conductor in a uniform space with a dielectric constant of ε and a magnetic permeability of μ. When it is affected by the incident field in any direction, it will generate a corresponding induced electromagnetic current and an outward scattered field. Assume that the incident electric field of the incident plane electromagnetic wave is E inc , the magnetic field is H inc , the scattered electric field generated by irradiating the surface of an ideal conductor is E sac , the scattered magnetic field is H sac , inside an ideal conductor, the internal electric field E in and the internal magnetic field H in All are zero; Step (2), assuming the simulation frequency band is set to [f1,f n ], according to the bisection interpolation theory, we know that there are three frequency points selected in the first round, namely f1, f n and f mid , f1 is the starting frequency point of the simulation, f n is the end frequency point of the simulation, f mid is the middle value between the start and end frequency points, f mid =(f1+f n ) / 2; Step (3), using the moment method to solve the electromagnetic field inside the conductor at the currently selected frequency point and constructing a matrix equation; Step (4) uses the accelerated GMRES solver to iteratively solve the matrix equation to obtain the current value at the simulation frequency point; specifically: 4-1 Calculate the current value at the current frequency point; specifically: 4-1-1 For a matrix equation (7), initialize I i =I0, I0 is initialized to 0, the number of iterations i = 0; 4-1-2 According to formula (7), we can know formula (11), and then get the residual r i ; r i =V-ZI i formula (11) 4-1-3 Determine the residual r i Is it less than the artificially set threshold k? If so, it is considered that the I i This is the current value you are looking for, otherwise go to step 4-1-4; 4-1-4 will be the residual r i The orthogonalized result is used as the basis vector of the initial Krylov subspace, and then Arnoldi iteration is performed to obtain an approximate solution vector I i+1 ; Then return to step 4-1-2 and update the number of iterations i = i + 1; 4-2 Calculate the current value at the third frequency point; specifically: 4-2-1 Calculate step (2) f1 and f according to step 4-1 n and f mid The current value of any two frequency points in the calculation is used to calculate the estimated current value of the third frequency point by interpolation algorithm. i Initialize to the estimated current value at the third frequency point, and initialize the number of iterations i = 0; 4-1-2 According to formula (11), the residual r is calculated i ; 4-1-3 Determine the residual r i If it is less than 1, then the estimated current value at the current frequency point is considered to be an effective value, otherwise I i Reset to 0 and repeat steps 4-1-2 to 4-1-4; Step (5), simulation accuracy evaluation: For step (2) f1, f n and f mid The current value result is used to make a convergence judgment. If the convergence condition is met, the simulation calculation is stopped. If the convergence condition is not met, the bisection interpolation theory is used to add new sampling frequency points for the next round, and steps (3)-(5) are repeated. Step (6): The electromagnetic simulation parameters can be obtained based on the current value obtained by the above calculation to realize electromagnetic simulation.

2. The method according to claim 1, characterized in that Step (3) is specifically: According to the electric field integral equation: In the formula represents the normal unit vector outside the conductor, η represents the wave impedance of the medium space, L represents the linear operator, J represents the current density, E inc represents the incident electric field; The RWG basis function is defined as being able to simulate the electromagnetic flow distribution on the surface of an object of arbitrary shape when there are common edges between two adjacent triangles; The definition of the RWG basis function is given as: where l n is the nth pair of adjacent triangles and The length of the common side, r is any point on the triangle, From the triangle Vertex The vector pointing to point r, From the triangle Vertex The vector pointing to point r, and Adjacent triangles and The area of ​​​​f n (r) represents the nth RWG basis function in the field region; By expanding the current in equation (1) using the RWG basis function: Among them I n is the nth current expansion coefficient to be solved, and N is the number of unknowns formed after the network is partitioned; Substituting formula (3) into formula (1) yields the discretized electric field integral equation: Where L(f n (r)) represents the RWG basis function f n (r) linear operator equation; The Galerkin method is used to select the RWG basis function as the test function to test the electric field integral equation (4): where f n (r) represents the nth RWG basis function in the field region; f m (r') represents the m-th RWG basis function at point r' in the source region; Since the RWG basis function is along the surface tangent direction, formula (5) can be rewritten as follows: Write formula (6) in the form of a matrix equation: ZI=V Formula (7) Where Z is the N×N moment method impedance matrix, V and I are N×1 column vectors, representing the excitation vector and current vector in the moment method, respectively; The elements of the mth row and nth column of the matrix Z and the mth row of the matrix V are expressed as follows: Where j represents the imaginary unit, f represents the frequency point; RWG basis function f m The source triangle pair where (r') is located is RWG basis function f n The field triangle where (r) is located is G(R) is the Green's function in the uniform space, R = |r-r'| represents the distance from any point r of the field basis function to any point r' of the source basis function; ▽ represents the gradient operator, ▽' represents the divergence operator, ds represents the surface vector element of the field region, and ds' represents the surface vector element of the source region.

3. The method according to claim 1, characterized in that The convergence criterion described in step (5) uses the error formula Where S(f) represents the theoretical value of frequency response, Represents the frequency response fitting value; if err < threshold, convergence is achieved, otherwise convergence is not achieved.

4. A computer-readable storage medium having a computer program stored thereon, which, when executed in a computer, causes the computer to execute the method according to any one of claims 1 to 3.

5. A computing device comprising a memory and a processor, wherein the memory stores executable code, and when the processor executes the executable code, the method according to any one of claims 1 to 3 is implemented.

Images