CN113887103A

CN113887103A - Integrated circuit full-wave electromagnetic simulation method and system based on different dielectric characteristics

Info

Publication number: CN113887103A
Application number: CN202111166781.8A
Authority: CN
Inventors: 王芬
Original assignee: Beijing Wisechip Simulation Technology Co Ltd
Current assignee: Beijing Wisechip Simulation Technology Co Ltd
Priority date: 2021-09-30
Filing date: 2021-09-30
Publication date: 2022-01-04
Anticipated expiration: 2041-09-30
Also published as: CN113887103B

Abstract

The application discloses an integrated circuit full-wave electromagnetic simulation method and system based on different dielectric characteristics, the method comprises the steps of firstly acquiring layer medium material information of the integrated circuit, determining the medium characteristics of the integrated circuit about conductivity, permeability and dielectric constant based on the layer medium material information, then, the medium type of the integrated circuit about loss and dispersion characteristics is judged based on the medium characteristics, a matrix equation corresponding to the medium type is established based on a mesh division result obtained after mesh division is carried out on a parallel flat plate field of the integrated circuit, a reference frequency point and a field solution thereof of the integrated circuit are calculated by an iteration method based on the matrix equation, and then, for a frequency point to be solved which is lower than the reference frequency point, and finally, obtaining the electromagnetic response of the full frequency band based on the field solutions of all the frequency points to be solved. The method can provide corresponding integrated circuit full-wave electromagnetic simulation schemes for different types of media.

Description

Integrated circuit full-wave electromagnetic simulation method and system based on different dielectric characteristics

Technical Field

The present application relates to the field of electromagnetic simulation technologies, and in particular, to a full-wave electromagnetic simulation method and system for an integrated circuit based on different dielectric characteristics.

Background

Very large scale integrated circuits have a pronounced multi-scale structure with dimensions in the centimeter (10) range^-2m) -nm class (10)^-9m), the scale range is up to 7 orders of magnitude. On the other hand, signals transmitted by integrated circuits often feature full-wave transmission, and the transmission frequency range covers from DC to higher harmonic frequencies, which is a problem in digital and mixed signal transmissionParticularly in integrated circuit applications. Therefore, the electromagnetic field analysis for the integrated circuit needs to be performed with a full-wave electromagnetic field analysis, and a full-wave electromagnetic field solver needs to be used to perform the full-wave electromagnetic field solution on the electromagnetic field of the integrated circuit.

However, investigation and test results show that, when a full-wave electromagnetic field is solved for the electromagnetic field problem of the integrated circuit, for the same integrated circuit model, when the test frequency is a high frequency above GHz, the solver can obtain an accurate field solution, and for the test frequency as low as MHz magnitude or lower than MHz magnitude, all solvers fail and cannot obtain an accurate field solution, and MHz and tens of MHz are just the working frequencies of many integrated circuits, so that the problem of solving the failure of the low-frequency electromagnetic field is urgently needed to be solved.

The reason for this is that, during the operation of the integrated circuit, the electromagnetic field forming the integrated circuit is formed by the superposition of the contributions of the conduction current and the displacement current, which are equivalent at high frequencies and therefore normally superposed; at low frequencies, however, the contribution of the conduction current is predominant, whereas the contribution of the displacement current is very low, the ratio of the contribution of the displacement current to the conduction current being even lower than the machine accuracy, e.g. at 10 for double-precision data storage^-16In order of magnitude, the ratio of the displacement current to the conduction current contribution at low frequencies will be comparable to or even lower than the machine accuracy, so that the ratio of the matrix elements representing the different contributions is of the order of magnitude comparable to or even lower than the machine accuracy. This fact causes that when the matrices representing different contributions are combined, the matrix elements representing the contributions of the displacement currents are completely covered by errors caused by machine precision during combination, so that the contributions of the displacement currents are negligible, but after the errors caused by machine precision are introduced, the contributions of the displacement currents are amplified by several orders of magnitude due to the errors, so that the contributions of the displacement currents are changed from negligible to resolvable, and the solution result is invalid, namely, the solved electromagnetic field is not an accurate field.

The existing method for solving the problem that the solver fails at low frequency is to combine the electromagnetic field solver based on constancy or quasi-constancy and the electromagnetic field solver based on high frequency to solve. When the test frequency is higher than a certain frequency, a high-frequency electromagnetic field solver is adopted for solving, and when the test frequency is lower than the certain frequency, a constant or quasi-constant electromagnetic field solver is adopted for solving, and then the calculation results of the two solvers are spliced.

However, this method is less accurate, firstly, because the constant or quasi-constant electromagnetic field solver involves a basic approximation, it is necessary to decouple the electric field E and the magnetic field H, forming a differential equation containing only the constant electric field or the constant magnetic field, which is correct only in the case of a strict dc field; secondly, how to set the specific frequency switched between the two solvers is unknown at present; finally, the field solved after being lower than a certain frequency is replaced by the constant field under direct current, so that the solved electromagnetic field is irrelevant to the frequency, the electromagnetic response curve of the integrated circuit of field splicing solved by the two solvers is obviously discontinuous, the electromagnetic response curve has obvious jump at the frequency point of switching of the two frequencies, and the response curve at a low frequency band is a straight line.

Therefore, in order to thoroughly solve the problem of calculation failure of the low-frequency field solution of the existing solver integrated circuit, it is necessary to accurately calculate and obtain the true solution of the full-wave maxwell equation set of the electric field E and the magnetic field H coupled from direct current to high frequency in the integrated circuit; moreover, how to accurately calculate and obtain the real solution under the condition that the dielectric characteristics of the integrated circuits are different is another problem which needs to be solved urgently.

Disclosure of Invention

Based on this, in order to solve the problem that the existing solver fails under the low-frequency condition of the integrated circuit under the condition of a lossless non-frequency-dispersion medium, and further obtain an electromagnetic field of the integrated circuit under the full frequency band including the low frequency, and accurately calculate and obtain the true solution of the integrated circuit with different medium characteristics, the application discloses the following technical scheme.

In one aspect, an integrated circuit full-wave electromagnetic simulation method based on different dielectric characteristics is provided, and includes:

acquiring layer medium material information of the integrated circuit, and determining medium characteristics of the integrated circuit, which relate to conductivity, magnetic permeability and dielectric constant, based on the layer medium material information;

determining a dielectric type of the integrated circuit with respect to loss and dispersion characteristics based on the dielectric characteristics;

establishing a matrix equation corresponding to the type of the medium based on a mesh generation result obtained after mesh generation is carried out on a parallel flat plate field of the integrated circuit;

calculating a reference frequency point and a field solution thereof of the integrated circuit by an iteration method based on the matrix equation, wherein the integrated circuit field solution corresponding to the reference frequency point is an accurate solution;

obtaining a field solution under a frequency point to be solved in a mode corresponding to the medium type, wherein when the medium type belongs to a lossy medium, the field solution of a reference frequency point is firstly blocked, the field solution of the frequency point to be solved is obtained according to the property of a finite element system matrix characteristic value and a blocking result, and when the medium type belongs to a dispersive medium, a proportionality coefficient between the field solution of the reference frequency point and the field solution of the low frequency point is firstly calculated, and the field solution under the frequency point to be solved is obtained according to the proportionality coefficient;

and obtaining the electromagnetic response of the full frequency band based on the field solutions of all the frequency points to be solved.

In a possible implementation manner, the calculating the reference frequency point of the integrated circuit by an iterative method based on the matrix equation includes:

step A1, setting an iteration frequency lower limit F_minIs a critical frequency point f₀And setting an upper iteration frequency limit F_max＝Factor×f₀Wherein Factor is the multiple of critical frequency point;

step A2, converting the current angular frequency omega_curr＝2πF_minSubstituting into a matrix equation corresponding to the integrated circuit dielectric loss type, and solving the matrix equation to obtain omega_currField solution at angular frequency E_curr(ii) a Wherein,

class of dielectric lossThe corresponding matrix equation when the model is a lossless medium is as follows: (K)₁-ω²K₂) E ═ b (ω), and when the dielectric loss type is a lossy medium, the corresponding matrix equation is: (K)₁-ω²K₂+jωK₃) E ═ b (ω), where ω is the electromagnetic angular frequency, j is the imaginary unit, E is the electric field, b (ω) is the external excitation source of the whole finite element system, K is₁Is the stiffness matrix of the entire finite element system, K₂Is a dielectric constant-dependent mass matrix, K, of the entire finite element system₃Is a mass matrix related to the electrical conductivity of the whole finite element system;

step A3, calculating the field solution E through a relative error calculation formula corresponding to the type of the dielectric loss of the integrated circuit_currRelative error res of; wherein, the calculation formula of the relative error res _ b under the lossless medium is as follows:

the formula for the relative error res _ a for lossy media is:

wherein, ω is_{curr_a}For angular frequency, omega, at lossy media_{curr_b}At angular frequency in a lossless medium, E_{curr_a}For field solutions in lossy media, E_{curr_b}Field decomposition under lossless medium;

step A4, when the relative error res is less than or equal to epsilon₁Then, a reference frequency point f is obtained_ref＝ω_currPer 2 pi and field solution E thereof_ref＝E_currAnd ending; at the relative error res>ε₁Jumping to step a 5;

step A5, mixing omega_curr＝π(F_min+F_max) Substituting the matrix equation corresponding to the integrated circuit dielectric loss type to obtain a new field solution;

step A6, calculating the relative error of the new field solution through the relative error calculation formula;

step A7, the relative error at the new field solution is less than ε₀When making F_max＝ω_currA/2 pi and jumping to step A5; the relative error at the new field solution is less than or equal to epsilon₁And is greater than or equal to epsilon₀Then, a reference frequency point f is obtained_ref＝ω_currAnd/2 pi and field solution thereof, and finishing; the relative error at the new field solution is greater than epsilon₁When making F_min＝ω_currA/2 pi and jumping to step A5; wherein epsilon₀Is a preset lower error threshold value epsilon₁Is a preset upper limit of the error threshold, epsilon₀<ε₁。

In a possible embodiment, when the medium type is a lossless non-dispersive medium, the E is determined according to a field solution at a reference frequency point_refCalculating a field solution e (f) for the frequency f to be solved by:

wherein f is<f_ref；

When the medium type is a lossless dispersive medium, obtaining a solution E (f) at a frequency point to be solved by the following formula: e (f) ═ kE_refWhere k is a proportionality coefficient and k is a real number, E_refIs a reference frequency point f_refThe following accurate field solutions, the proportionality coefficient k being:

wherein,

is a non-zero vector.

In one possible embodiment, when the media type is lossy, non-dispersive media or lossy, dispersive media:

before the matrix equation corresponding to the medium type is established, dividing materials for forming the multilayer integrated circuit into non-conductive media and conductive media, numbering grid nodes based on the split parallel flat plate field according to the medium areas, enabling the grid nodes to be continuously arranged according to the medium areas where the grid nodes are located, and taking the results of numbering and arranging as the grid split result; and,

before calculating the reference frequency point and the field solution of the integrated circuit, dividing the finite element unknowns according to the number of the grid nodes in the non-conductive medium area and the number result of the grid nodes, and then partitioning the stiffness matrix of the equation set according to the sorting result of the finite element unknowns, wherein the stiffness matrix after partitioning is the stiffness matrix after partitioning

Wherein n is a finite element unknown quantity in the non-conductive medium region, c is a finite element unknown quantity in the conductive medium region, and the finite element unknown quantity is an electric field of the basic unit edge or surface element; k (omega) is the left-end matrix of the finite element system matrix, K_nn(ω) is a submatrix associated with a region of a non-conductive medium, K_nc(ω) is a submatrix of non-conductive medium regions associated with conductive medium regions, K_cn(ω) is a submatrix of regions of conductive medium associated with regions of non-conductive medium, K_cc(ω) is a submatrix associated with a region of conductive medium; and,

the blocking of the field solution of the reference frequency point comprises the following steps: partitioning the field solution of the reference frequency point according to the partitioning result of the stiffness matrix of the equation set to obtain

Wherein E is_n,refFor the part of the field solution at the reference frequency lying within the non-conductive medium region, E_c,refIs the part of the field solution under the reference frequency point in the conductive medium area.

In one possible embodiment, when the media type is lossy, frequency-dispersive media: the acquisition of the field solution under the frequency point to be solved comprises the following steps:

splitting the field solution of the reference frequency point into a real part and an imaginary part, and setting the imaginary part in the conductive medium area to zero to obtain the split field solution of the reference frequency point:

wherein re (·) is a real part, im (·) is an imaginary part, j is an imaginary unit, and then the field solution of the frequency to be solved is calculated according to the relationship between the field solution of the frequency to be solved and the split reference frequency point field solution by the following formula:

when the media type is lossy dispersive media: the acquisition of the field solution under the frequency point to be solved comprises the following steps:

wherein E is_n,ref,reFor field solutions at the reference frequency point lying in the real part, E, of the non-conductive medium region_n,ref,imFor field solutions at the reference frequency point in the imaginary part, E, of the non-conducting medium region_c,ref,reThe field solution of the frequency to be solved is located in the real part of the conductive medium area under the reference frequency point, and then is calculated according to the relationship between the field solution of the frequency to be solved and the split reference frequency point field solution by the following formula: e (f) zE_refWherein z is a proportionality coefficient and z is a complex number, the proportionality coefficient z being equal to z_re+jz_imWherein z is_reIs the real part of the proportionality coefficient, z_imIs the imaginary part of the scaling factor;

z is_reAnd z_imIs given by the following formula:

wherein,

matrix coefficients that are the real part of the scaling coefficient z,

matrix coefficients being the imaginary part of the scaling coefficient z,

is E_n,ref,imThe transpose matrix of (a) is,

is E_n,ref,reTransposed matrix of b₀Is a constant vector independent of frequency;

the above-mentioned

And

obtained by the following formula:

wherein, K_2,nnIs a matrix K₂Submatrix, K, associated with areas of non-conducting medium_3,ccIs a matrix K₃With a sub-matrix associated with the region of conductive medium.

In another aspect, an integrated circuit full-wave electromagnetic simulation system based on different dielectric characteristics is further provided, including:

the medium information acquisition module is used for acquiring layer medium material information of the integrated circuit and determining medium characteristics of the integrated circuit, which relate to conductivity, permeability and dielectric constant, based on the layer medium material information;

the medium type judging module is used for judging the medium type of the integrated circuit about loss and dispersion characteristics based on the medium characteristics;

the system comprises an equation set building module, a matrix equation generating module and a data processing module, wherein the equation set building module is used for building a matrix equation corresponding to the type of the medium based on a mesh division result obtained by carrying out mesh division on a parallel flat plate field of the integrated circuit;

a reference frequency point obtaining module, configured to calculate a reference frequency point and a field solution thereof of the integrated circuit through an iterative method based on the matrix equation, where the integrated circuit field solution corresponding to the reference frequency point is an accurate solution;

the to-be-solved field solution calculation module is used for acquiring field solutions under frequency points to be solved in a mode corresponding to the medium type, wherein when the medium type belongs to a lossy medium, the field solutions of the reference frequency points are firstly partitioned, the field solutions under the frequency points to be solved are acquired according to the property of the matrix characteristic value of the finite element system and the partitioning result, and when the medium type belongs to a dispersive medium, the proportionality coefficient between the field solutions of the reference frequency points and the low-frequency field solutions is firstly calculated, and the field solutions under the frequency points to be solved are acquired according to the proportionality coefficient;

and the electromagnetic response acquisition module is used for acquiring the electromagnetic response of the full frequency band based on the field solutions of all the frequency points to be solved.

In a possible implementation manner, the reference frequency point obtaining module calculates the reference frequency point of the integrated circuit by the following steps:

when the medium loss type is a lossless medium, the corresponding matrix equation is as follows: (K)₁-ω²K₂) E ═ b (ω), and when the dielectric loss type is a lossy medium, the corresponding matrix equation is: (K)₁-ω²K₂+jωK₃) E ═ b (ω), where ω is the electromagnetic angular frequency, j is the imaginary unit, E is the electric field, b (ω) is the external excitation source of the whole finite element system, K is₁Is the stiffness matrix of the entire finite element system, K₂Is a dielectric constant-dependent mass matrix, K, of the entire finite element system₃Is a mass matrix related to the electrical conductivity of the whole finite element system;

the formula for the relative error res _ a for lossy media is:

In one possible embodiment, inWhen the medium type is a lossless non-frequency dispersion medium, the to-be-solved field solution calculation module performs field solution E under a reference frequency point_refCalculating a field solution e (f) for the frequency f to be solved by:

wherein f is<f_ref；

When the medium type is a lossless dispersive medium, the to-be-solved field solution calculation module obtains a solution E (f) at a to-be-solved frequency point by the following formula: e (f) ═ kE_refWhere k is a proportionality coefficient and k is a real number, E_refIs a reference frequency point f_refThe following accurate field solutions, the proportionality coefficient k being:

wherein,

is a non-zero vector.

before the equation structure modeling block establishes a matrix equation corresponding to the type of the medium, dividing a material for forming a multilayer integrated circuit into a non-conductive medium and a conductive medium, numbering grid nodes based on a subdivided parallel flat plate field according to the medium area, continuously arranging the grid nodes according to the medium area where the grid nodes are located, and taking the numbering and arranging results as the mesh subdivision results; and,

before the reference frequency point acquisition module calculates the reference frequency point and the field solution of the integrated circuit, the finite element unknowns are divided according to the number of the grid nodes in the non-conductive medium area and the number result of the grid nodes, and then the stiffness matrix of the equation set is partitioned according to the sorting result of the finite element unknowns, wherein the stiffness matrix after partitioning is the stiffness matrix after partitioning

the to-be-solved field solution calculation module divides the field solution of the reference frequency point into blocks by the following steps:

partitioning the field solution of the reference frequency point according to the partitioning result of the stiffness matrix of the equation set to obtain

In one possible embodiment, when the media type is lossy, frequency-dispersive media: the to-be-solved field solution calculation module acquires a field solution under a to-be-solved frequency point through the following steps:

splitting the field solution of the reference frequency point into a real part and an imaginary part, and zeroing the imaginary part in the conductive medium area to obtain the split field solution of the reference frequency point:

when the media type is lossy dispersive media: the to-be-solved field solution calculation module acquires a field solution under a to-be-solved frequency point through the following steps:

z is_reAnd z_imIs given by the following formula:

wherein,

matrix coefficients that are the real part of the scaling coefficient z,

matrix coefficients being the imaginary part of the scaling coefficient z,

is E_n,ref,imThe transpose matrix of (a) is,

the above-mentioned

And

obtained by the following formula:

According to the integrated circuit full-wave electromagnetic simulation method and system based on different dielectric characteristics, the types of the media can be divided based on different dielectric characteristics, and corresponding integrated circuit full-wave electromagnetic simulation schemes are provided for the different types of media.

In addition, the method can calculate the reference frequency point and the solution under the reference frequency point by adopting an iteration method, and reversely deduces the solution under the low frequency based on the reliable reference frequency point which can obtain an accurate solution, so as to obtain the electromagnetic field of the integrated circuit under the low frequency, thereby solving the failure problem of the existing solver under the low frequency condition aiming at the integrated circuit, realizing the accurate solution of the field solution of the low frequency point, simultaneously solving the continuity problem of low frequency and high frequency response, leading the simulation result to be more accurate, and avoiding the problem of discontinuous curves of response splicing of two solvers caused by the difference of the solution results of the frequency points where the high frequency band and the low frequency band intersect when different solvers are respectively adopted aiming at the high frequency and the low frequency.

In addition, the method and the device utilize the characteristic that the real part of the electromagnetic response of the integrated circuit under the lossy medium is irrelevant to the frequency of the field at the low-frequency point and the imaginary part of the electromagnetic response is in inverse proportion to the frequency, based on the characteristic, the field solution under the low frequency is reversely deduced through the obtained accurate field solution under the reliable reference frequency point under the lossy medium, so that the electromagnetic response of the integrated circuit under the low frequency is obtained, meanwhile, the characteristic that the expansion characteristic exists between the field solution under the low frequency to be solved under the lossy medium and the accurate field solution under the reference frequency point is utilized, based on the characteristic, the field solution under the low frequency is reversely deduced through the obtained accurate field solution under the reliable reference frequency point under the lossy medium, and further the electromagnetic response of the integrated circuit under the low frequency is obtained.

Drawings

The embodiments described below with reference to the drawings are exemplary and intended to be used for explaining and illustrating the present application and should not be construed as limiting the scope of the present application.

FIG. 1 is a schematic flow diagram of an embodiment of a method for full-wave electromagnetic simulation of an integrated circuit based on different dielectric properties disclosed in the present application.

FIG. 2 is a block diagram of an embodiment of an integrated circuit full-wave electromagnetic simulation system based on different dielectric characteristics.

Detailed Description

In order to make the implementation objects, technical solutions and advantages of the present application clearer, the technical solutions in the embodiments of the present application will be described in more detail below with reference to the drawings in the embodiments of the present application.

An embodiment of the integrated circuit full-wave electromagnetic simulation method based on different dielectric characteristics disclosed in the present application is described in detail below with reference to fig. 1.

As shown in fig. 1, the method disclosed in the present embodiment includes the following steps 100 to 600.

Step 100, layer dielectric material information of the integrated circuit is obtained, and dielectric properties of the integrated circuit, which relate to conductivity, permeability and permittivity, are determined based on the layer dielectric material information.

Specifically, material information of each layer of medium defined by the multilayer integrated circuit model is obtained from an object for storing the integrated circuit model, the material information comprises the thickness of each layer of medium layer and the thickness of copper-coated layer, the electric conductivity, the magnetic conductivity and the dielectric constant of each layer of medium layer at different temperatures are obtained aiming at the electric-thermal combined simulation of the integrated circuit, and the electric conductivity, the magnetic conductivity and the dielectric constant of each layer of medium layer at different frequencies are obtained aiming at the full-wave electromagnetic simulation of the integrated circuit.

Step 200, judging the medium type of the integrated circuit about loss and dispersion characteristics based on the medium characteristics.

Specifically, the conductivity of the dielectric layer obtained in step 100 is determined, and if the conductivity of the dielectric layer is greater than 0, the type of the dielectric of the integrated circuit is determined to be a lossy dielectric, but if the conductivities of all the dielectric layers are 0, that is, the dielectric layer used by the integrated circuit is considered to be an ideal insulating dielectric layer, the type of the dielectric of the integrated circuit is determined to be a lossless dielectric. In addition, the change conditions of the conductivity, the permeability and the permittivity of the dielectric layer obtained in step 100 at different frequencies are also determined, if the conductivity, the permeability and the permittivity are not changed, it is determined that the dielectric layer used by the integrated circuit does not have the frequency dispersion characteristic, but if at least one of the conductivity, the permeability and the permittivity is changed, it is determined that the dielectric layer used by the integrated circuit has the frequency dispersion characteristic.

Therefore, the medium types can be divided according to the existence of loss and the existence of frequency dispersion, so that the medium types are four in total, namely, firstly, the medium without loss and the frequency dispersion, secondly, the medium with loss and the frequency dispersion, thirdly, the medium without loss and the frequency dispersion and fourthly, the medium with loss and the frequency dispersion.

In one embodiment, after step 200, the following steps B1 and B2 may be performed before step 300 is performed.

And step B1, obtaining model information of the very large scale integrated circuit and modeling the integrated circuit according to the model information. The model information comprises layer information, layout information of each layer, via hole information and netlist information of the integrated circuit.

When performing electromagnetic simulation of the very large scale integrated circuit, the chip design file may be read first to obtain the model information of the very large scale integrated circuit, that is, the layer information of the integrated circuit, the layout information of each layer of the integrated circuit, the via information of the integrated circuit, and the netlist information of the integrated circuit.

Wherein the layer information of the integrated circuit may include: number of layers, layer thickness, dielectric material information of each dielectric layer, material and thickness of each conductive layer, etc.

Regarding layout information of each layer of the integrated circuit, the layout information of each layer of the integrated circuit is generally described by adopting polygons, the polygons are generally formed by connecting a series of vertexes, the sequence of the vertexes of the polygons is positive anticlockwise, a conductive power supply layer is filled in the polygons, the polygons are negative clockwise, the polygons are hollowed polygons, and insulating media are arranged in the polygons. The routing of the integrated circuit is converted into polygonal information in the layout of the integrated circuit.

The via information of the integrated circuit may include: via position, via size, and layer number of the integrated circuit to which the upper and lower vertices of the via are connected. The through holes of the integrated circuits are through holes for connecting the two layers of integrated circuits, and the two layers of integrated circuits are electrically connected at a proper position by introducing the through holes at the position.

With regard to the netlist information of the integrated circuit, the connection relationship of the external circuit of the integrated circuit, that is, topology information, and information of components, power supplies and the like constituting the circuit can be acquired through the netlist information. The circuit information and the network information designed on each layer of the layout of the integrated circuit can be obtained through the netlist information, wherein the circuit information provides a series of node information and connection relations thereof which are positioned on different layers and different positions, the network information and the circuit information are in an inclusion relation, and one network can contain a plurality of circuits. The above node information, circuit information and network information are given by unique names.

After the model information is obtained based on the chip design file, the integrated circuit can be modeled by using the model information of the integrated circuit, so as to obtain an integrated circuit model. The modeling process may include: converting layout information of the integrated circuit into discrete point information and polygon information containing layer information formed by the discrete points; converting the via hole information of the integrated circuit into a polygon prism (such as a hexagonal prism and a twelve-prism) connecting two layers; and converting the netlist information of the integrated circuit into external circuit information of the integrated circuit and topological structure information of the layout node.

And step B2, performing mesh generation on the parallel flat plate field of the integrated circuit based on the model obtained by modeling.

In the case of an alternating electromagnetic field, electromagnetic waves in the integrated circuit propagate through a medium between metal layers of different layers, and a medium region between the metal layers of different layers is called a parallel flat plate field, so that in the simulation calculation of the alternating electromagnetic field of the integrated circuit, the parallel flat plate field needs to be identified from the model obtained by modeling in the step B1 and then is subjected to mesh subdivision.

After the modeling of the multilayer integrated circuit is completed, the layout polygons of the multilayer integrated circuit can be simplified, and then mesh generation can be performed on the modeling model after layout simplification. The simplification is to make the subsequent mesh generation process more accurate and rapid, specifically, because the layout information in the chip design file is the design file of the actual physical model, the original file may be a graphic file, after being converted into a layout polygon, a large amount of redundant node information may exist, when the electromagnetic field numerical calculation method is adopted to perform electromagnetic simulation on the integrated circuit, the calculated field needs to be dispersed based on the structure of the integrated circuit, namely mesh generation, and then a discrete equation set is established based on the generated mesh to solve.

If the simulation is directly performed on layout polygon information in a chip design file, nodes of each input polygon are regarded as important nodes, so that all input polygon nodes are inserted into split mesh nodes, which may generate very dense meshes with poor quality, so that a large amount of computer resources are wasted while inaccurate calculation is performed. Therefore, before mesh subdivision is carried out on the layout polygons of the multilayer integrated circuit, each layout polygon of each layer can be simplified without losing precision, so that compared with the original input polygons, the shapes of the simplified polygons do not change within the precision control range, and most importantly, the bad consequences that the connection relation of the polygons changes due to the change of the polygon shapes does not occur, and the topological structure of the integrated circuit layout is changed are avoided. After polygon reduction, the vertices of the polygons reaching the minimum number form polygons of the same shape as the original input polygon.

After the layout is simplified, the alignment of the polygons of the layout can be performed before the mesh generation, and then the mesh generation can be performed on the parallel flat plate field after the polygons are aligned. Specifically, if two layers of integrated circuit layouts designed have a layout polygon that is the same and completely aligned, a parallel flat field with the shape of the layout polygon will be formed between the two layers of layouts, however, since actually provided layout polygon information may generate some errors due to the processes of file format conversion and the like, the errors may cause that the finally introduced layout polygons of the two layers are not completely the same or not completely aligned, and compared with the polygon size, there may be an error of 1% or 0.1%. If the parallel flat plate field is directly identified based on the input polygon information, a plurality of parallel flat plate field fragments are generated, the fragment size is very small, the introduction of the fragment may also cause very dense mesh subdivision with poor quality, the result is wrong while the computing resource is wasted, and therefore, the layout polygons of the multilayer integrated circuit can be aligned without losing precision so as to obtain the parallel flat plate field with higher accuracy and compactness.

After the alignment of the multilayer integrated circuit layout is realized, the parallel flat plate field of the multilayer integrated circuit layout can be identified based on the layer information of the polygon of the integrated circuit layout.

And 300, establishing a matrix equation corresponding to the type of the medium based on a mesh generation result obtained after mesh generation is carried out on the parallel flat plate field of the integrated circuit.

In one embodiment, step 300 specifically includes the following steps 310 to 350.

And step 310, establishing an electromagnetic field wave equation based on the Maxwell equation set. Specifically, step 310 first establishes a wave equation based on the electric field E using maxwell's equations to obtain an electromagnetic field wave equation in the following formula (1):

in equation (1),. DELTA.is the rotation operator,. mu._rIn order to obtain a relative magnetic permeability of the medium,

is the electric field vector, ω is the electromagnetic angular frequency (in rad/s), c0 is the wave velocity of the electromagnetic wave in vacuum, c0 is 3 × 10⁸m/s，ε_rIs the relative dielectric constant of the medium, j is an imaginary unit, j²＝-1，μ₀Is the permeability of a vacuum medium, mu₀＝4π×10^-7H/m, σ is the conductivity of the medium (in S/m),

current density (in A/m) for applied excitation²)。

And 320, acquiring a homogeneous equation corresponding to the electromagnetic field wave equation to obtain a functional of the homogeneous equation. Specifically, in step 320, a homogeneous equation corresponding to the electromagnetic field fluctuation equation is obtained by a variational principle, and the functional is expressed by the following formula (2):

in the formula (2), V is an electromagnetic field resolving area, S is a plane surrounded by the electromagnetic field resolving area, and n is a normal vector in which any point of the plane S points outward.

And step 330, when the size of the electromagnetic field solving area reaches a set threshold value, setting the part of the functional, which is related to the electromagnetic wave and is at the boundary of the area, as 0. Specifically, when the electromagnetic field solution area is large enough in step 330, so that the electromagnetic wave is attenuated to approximately 0 at the boundary of the area, the functional function can be simplified as the following formula (3):

in the equation (3), the judgment method of whether the electromagnetic field solution area is large enough may be to obtain a ratio between a minimum distance from a boundary of the solution area to a source (i.e., a multilayer integrated circuit board) generating the electromagnetic wave and a wavelength of the electromagnetic wave, compare the ratio with a preset multiple, and determine that the electromagnetic field solution area is "large enough" if the ratio exceeds the preset multiple, for example, the preset multiple is 10 times, and determine that the solution area is large enough if the ratio is greater than 10.

And 340, dispersing the electromagnetic field solving area based on the basic geometric units formed by mesh subdivision, and interpolating the electric field of any point in each discrete unit by using the electric field of the edge or surface element of each discrete unit to obtain a discrete form functional of the functional shown in the formula (3). Specifically, in step 340, after the functional is simplified to obtain equation (3), the electromagnetic field solution area is discretized by using a basic unit which is small enough, where the basic unit may be a tetrahedron, a triangular prism, a hexahedron, or the like, and the electric field of any point in each discrete unit is interpolated by using an interpolation basis function and an electric field of an edge or a surface element, where the interpolation function used may specifically be equation (4) below:

in the formula (4), the reaction mixture is,

is the electric field of any point in the basic unit e, M is the number of interpolation basis functions,

the ith interpolation basis function of the basic unit E, E^eAn electric field vector formed by an electric field at an edge of the basic cell e,

the electric field value of the corresponding edge or surface element of the ith basis function of the basic unit e,

is a rib of a basic unitM interpolation basis functions on an edge or bin

Of size M × 1, { E^eIs M interpolation basis functions on the edge or surface element of the basic unit

Corresponding electric field value

The size of the matrix form of (1) is M × 1, and T represents a transpose of the matrix.

The determination method of whether the basic unit is small enough needs to be determined according to the following two conditions, and if the following two conditions are satisfied simultaneously, the size of the basic unit is considered small enough:

1. judging that the size of the nearby grid is not larger than the characteristic size of the integrated circuit according to the size relation between the local size of the basic unit and the characteristic size of the integrated circuit;

2. the relationship between the maximum size of the basic unit and the minimum wavelength of the electromagnetic wave of the integrated circuit to be simulated, wherein the minimum wavelength of the electromagnetic wave of the integrated circuit is not less than a preset multiple, such as 10 times, of the maximum size of the basic unit.

Substituting the interpolation function into the simplified equation (3) to obtain a discrete functional expressed by the following equation (5):

in the formula (5), the reaction mixture is,

wherein,

is the relative permeability of the medium in the region of the basic element e,

is the relative permittivity of the medium in the region of the elementary cell e,

p-th interpolation basis function of basic unit e, V^eIs an integral of the basic unit e, and L is the number of the basic units e in which the entire electromagnetic field resolving area is dispersed.

Is a stiffness matrix of the basic cell e,

is a dielectric constant dependent mass matrix of the medium of the elementary cell e,

is the conductivity-dependent mass matrix of the medium of the elementary cell e.

When an external stimulus is present, the formula (5) becomes the following formula (9):

in the formula (9), b^eIs an external excitation source of the basic unit e.

In the formula (10), the compound represented by the formula (10),

an external excitation current source, V, for the basic cell e^eIs an integral of the basic cell e.

And 350, taking a partial derivative of the discrete functional and setting the partial derivative to be 0 to obtain a matrix equation of the vector finite element under the lossless frequency dispersion-free medium. Specifically, in step 350, according to the energy minimization principle, the solution corresponding to the electromagnetic field wave equation (1) is the electric field E corresponding to the extremum value calculated by the functional expressed by the above equation (5), so that the following equation (11) can be obtained by taking the partial derivative of the equation (5) and making the partial derivative be 0:

(K₁-ω²K₂+jωK₃)E＝0 (11)；

in the formula (11), K₁Is the stiffness matrix of the entire finite element system, K₂Is a dielectric constant-dependent mass matrix, K, of the entire finite element system₃Is the conductivity-related mass matrix of the entire finite element system.

In the presence of an applied stimulus, a matrix equation is thus constructed as shown in equation (12 a):

(K₁-ω²K₂+jωK₃)E＝-b(ω) (12a)；

wherein b (omega) is an external excitation source of the whole finite element system and is a function of the angular frequency omega of the electromagnetic wave.

For the third lossy non-dispersive medium and the fourth lossy dispersive medium, both of which belong to lossy media, under the condition of lossy media, the conductivity of the dielectric layer of the integrated circuit is a non-zero finite value; in the dispersion type, the conductivity, dielectric constant and permeability of the dielectric layer and the metal layer of the integrated circuit belong to a dispersion-free medium, and the dielectric layer and the metal layer of the integrated circuit under the dispersion-free medium do not change along with the frequency or change along with the frequency is extremely small to be ignored, namely K is₁、K₂And K₃The dielectric constant and the magnetic permeability of the dielectric layer of the integrated circuit belong to a frequency dispersion medium, and the dielectric constant and the magnetic permeability of the dielectric layer of the integrated circuit under the condition of the frequency dispersion medium are all changed along with the change of the frequencyThe matrix K in the formula (12a) changes due to the change of the rate₁、K₂And K₃The element(s) of (c) will change due to the change in frequency. However, for both (c) and (c), the matrix equations are (12a), where the stiffness matrix corresponding to the media types of (c) and (c) is K in equation (12a)₁The corresponding mass matrix is K in the formula (12a)₂And K₃。

For a lossless non-frequency dispersion medium and a lossless frequency dispersion medium, the lossless non-frequency dispersion medium and the lossless frequency dispersion medium belong to lossless media, under the condition of the lossless media, a metal layer of an integrated circuit is an ideal conductor, and the conductivity of a dielectric layer is 0; in the dispersion type, the influence of the dispersion type is the same as that of the dispersion type in the third and fourth, and the description thereof is omitted. Thus, for (r) and (r) media whose media type is lossless, the stiffness matrix corresponding to the media type is K in equation (12a)₁The corresponding mass matrix is K in the formula (12a)₂However, since the conductivity is 0 at this time, K is₃0. Thus, the matrix equation corresponding to (r) and (c) is actually the following equation (12 b):

(K₁-ω²K₂)E＝-b(ω) (12b)。

and step 400, calculating the reference frequency point and the field solution of the integrated circuit by an iteration method based on the matrix equation. The reference frequency point means that the field solution of the integrated circuit calculated at the reference frequency point is definitely accurate, the direct current eigenmode of the electromagnetic wave at the reference frequency point plays a main role, and the contribution of the high-order eigenmode can be ignored.

In one embodiment, step 400 may include step 410 and step 420.

And step 410, calculating critical frequency points of the full-wave electromagnetic analysis of the integrated circuit based on the layout characteristic dimension of the integrated circuit and the machine precision of the simulation operation equipment. The critical frequency point is a critical frequency point, when a matrix equation for electromagnetic field simulation of the integrated circuit is solved from high frequency to low frequency in a frequency domain, the error of a solving result changes from small (result credibility) to large, when the error is large to a certain degree, the solving result is not credible, and the critical frequency from credible solving result to incredible solving result is the critical frequency point.

The characteristic dimension of the layout refers to the maximum dimension and the minimum dimension of the layout, the maximum dimension refers to the overall dimension of the whole layout plane, the minimum dimension refers to the dimension of the minimum functional unit of the layout, and the minimum dimension may be the diameter of the minimum via hole in the layout, the width of the thinnest routing, the width of the minimum gap between the routing, or the minimum layer thickness of the multilayer integrated circuit layout.

In one embodiment, step 410 may include step 411 and step 412.

And 411, obtaining magnitude ratios related to sizes among different matrix elements based on the range of the layout characteristic size.

In the most advanced very large scale integrated circuits at present, the feature sizes of the layouts of different positions of the integrated circuit have dimensions in the centimeter scale (10)^-2m) -nm class (10)^-9m), if tetrahedron is used to disperse the calculation field of the multilayer VLSI, the smallest dimension of the dispersed tetrahedron grid is nanometer.

In addition, in the matrix expressions in expressions (6) to (8), all the basic units

Is proportional to 1/l because of the interpolation function of the basic unit

Normalization has been performed, i is the size of the elementary cell obtained by mesh division, and the volume of the elementary cell e and i³Is in direct proportion, and

and

are all constants independent of layout feature size, so K in equation (6)₁Of the order of O (l), O (-) represents the acquisition order,o (l) represents the magnitude of l, K in formula (7)₂Of the order of O (c 0)^-2l³) K in the formula (8)₃Of the order of O (mu)₀σl³) Thus, the matrix K in the formula (13) can be obtained₁、K₂Norm ratio of (2):

and obtaining a matrix K in formula (14)₁、K₃Norm ratio of (2):

the formulae (13) and (14) can embody the matrix K₁、K₂And K₃There is inevitably a difference in order of magnitude between the elements (a) and (b), and it is also known from equations (13) and (14) that the discrete cell size (l) is the only change matrix (K)₁、K₂And K₃And the smaller the size l, the smaller the matrix K₁、K₂And K₃The larger the element contrast difference, and the smaller the difference. And due to the fact that

Therefore, consider the matrix K₁、K₂And K₃When the element contrast difference is a factor, only K is generally considered₁Norm and K of₂That is, if K is solved₁Norm and K of₂The problem of too large a norm of (A) naturally solves K₁Norm and K of₃Too large a ratio of norms of (a).

It can be understood that the layout characteristic size determines the distribution of the mesh subdivision size, for example, the layout is a multi-scale structure, and the characteristic size is from centimeter level of the maximum size to nanometer level of the minimum size, so that the maximum size of the subdivided tetrahedral unit may be centimeter level and is distributed at the position without the small-size layout; the minimum size is nano-scale and is distributed at the position of the small-size layout. Step B2 is therefore to subdivide the grid based on the layout information of the integrated circuit, the minimum feature size of the integrated circuit determining the minimum size of the grid elements, and the magnitude difference between the finite element system matrices in this step is based on the minimum size of the grid elements.

Suppose that in the integrated circuit currently implementing the method, the dimension range of the characteristic dimension of the layout of the integrated circuit at different positions is centimeter level (10)^-2m) -nm class (10)^-9m), if tetrahedron is used as basic unit to mesh the calculation field of multilayer VLSI, i.e. to disperse, then no matter for lossless media (r) and (r) or lossy media (r) and (r), the order of c0 is 10⁸L is of the order of 10^-9Thus K₂Norm and K of₁The ratio of the norm of (a) will also be as low as 10^-34Of the order of magnitude, i.e. a magnitude ratio of 10³⁴(ii) a For the lossy media (c) and (c) alone, at MHz, due to μ₀Of the order of 10^-7Assuming a value of 10 for σ³Then K is₁Norm and K of₃The ratio of the norm of (a) will also be as low as 10^-22Of the order of magnitude, i.e. a magnitude ratio of 10²²。

Step 412, obtaining the machine precision adopted in the simulation operation, and calculating the critical frequency point of the integrated circuit based on the machine precision and the magnitude ratio.

As can be seen from the stiffness matrix equations shown in equations (12a) and (12b) and the preceding discussion, all solvers fail at low frequencies, resulting in failure of the electromagnetic field simulation solver for integrated circuits due to the limited machine accuracy of the computer due to the matrix K in equations (12a) and (12b) formed using vector finite element calculations₁、K₂And K₃Is different by an order of magnitude which is inevitable when the frequency is low enough to make the frequency-dependent matrix ω in equation (12a) or equation (12b)²K₂Or j ω K₃Can lead to failure of the solver when the contribution of (a) is lost due to limited machine accuracy, since the left hand table of the matrix equation at this timeHas an expression approximately equal to K₁When this occurs, the solutions of equations (12a) and (12b) solved by the solver are completely wrong, because K is now the case₁Is a singular matrix.

Assuming that double-precision type data is adopted for calculation at present, the precision of a machine adopted in simulation operation is 10^-16When using double precision data types to perform omega²K₂And K₁In the subtraction of (2), if the matrix K is₁And ω²K₂The phase difference is greater than 10¹⁶The simulation computing device (e.g. computer) will directly use the matrix ω²K₂The value is zero, and the solution of equations (12a) and (12b) of the electromagnetic field simulation of the integrated circuit must fail. Even if the frequency is MHz, omega²K₂Is also compared with K₁To be smaller by 34-2 x 6-22 orders of magnitude, ω K₃Ratio K₁22-6 to 16 orders of magnitude smaller, so when ω is performed²K₂And K₁Is subtracted and ω K₃And K₁In the subtraction, the simulation computing device (e.g. computer) will directly apply the matrix ω²K₂And ω K₃Considered as zero.

Therefore, the cutoff frequency (namely, the critical frequency point) which fails in the simulation of the integrated circuit needs to be calculated based on the machine precision and the magnitude ratio, namely, the frequency which meets the following formula is the critical frequency point:

wherein f is the frequency point to be solved which is not higher than the critical frequency point. From this formula, one can obtain:

therefore, the formula for calculating the critical frequency point is as follows (15):

thus, a critical frequency point calculation formula is obtained, wherein f₀The critical frequency points are a machine precision magnitude adopted in simulation operation, c0 is the wave velocity of electromagnetic waves in vacuum, and l is the size of a basic unit obtained by mesh division.

Assuming that the machine precision magnitude a is 16, the formula for determining whether the critical frequency point condition is satisfied is specifically: II | K₁‖/‖(2πf)²K₂‖>10¹⁶The critical frequency point is

Due to the ratio II K₁‖/‖K₂Is difficult to be directly and accurately calculated, so the O (·) function is adopted to obtain the ratio II K₁‖/‖K₂Of the order of magnitude II, by O (c 0)²/l²) Instead of | K₁‖/‖K₂II this ratio is followed, but because of O (c 0)²/l²) Obtained is only c0²/l²Rather than an exact value, it is c0 that has an exact value²/l²Thus making it possible to

In one embodiment, l may be the smallest possible dimension, i.e., l ═ l_min，l_minIs the smallest dimension in the range of values for l. The value of l can be various, but for the most advanced super large scale integrated circuit, the minimum characteristic dimension of the layer structure and the layout reaches the nanometer level (10)^-9m), the size of the discrete tetrahedrons is also in the nanometer range, while if a grid size in the nanometer range can be achieved, a grid size higher than that of the nanometer process can be achieved, so l can be taken_min＝10^-9m, i.e. l 10^- ⁹m, critical frequency point f at this moment₀160MHz, that is, the results obtained when solving the integrated circuit electromagnetic field simulation matrix equations at frequencies below 160MHz must be inaccurate.

Steps 411 and 412 can calculate critical frequency points of full-wave analysis based on the feature size of the integrated circuit and the solution accuracy of the computer, so that the electromagnetic field simulation field solution at which frequency is obtained is definitely inaccurate. After the critical frequency point is obtained, the reference frequency point can be calculated according to the critical frequency point.

And step 420, calculating the reference frequency point and the field solution thereof by an iterative calculation method based on the critical frequency point.

In one embodiment, step 420 includes steps 421 through 427.

Step 421, setting the lower limit F of the iteration frequency_minIs a critical frequency point f₀And setting an upper limit initial value F of the iteration frequency_max＝Factor×f₀Wherein Factor is the multiple of the critical frequency point. Factor>Factor may be set to 10, 1. Specifically, for the case that the medium of the integrated circuit is a lossy medium or a lossless medium, setting respective lower limit initial value and upper limit initial value of iteration frequency, wherein F_{min_a}And F_{max_a}Respectively setting the lower limit initial value and the upper limit initial value of the iteration frequency under the lossy media tri and tetra, and setting F_{min_a}＝F_min，F_{max_a}＝F_max；F_{min_b}And F_{max_b}Respectively setting an iteration frequency lower limit initial value and an iteration frequency upper limit initial value under the lossless mediums of the first step and the second step, and setting F_{min_b}＝F_min，F_{max_b}＝F_max。

Step 422, the current angular frequency ω is calculated_curr＝2πF_minSubstituting the matrix equations of the formula (12a) and the formula (12b), and solving the matrix equations to obtain omega_currField solution at angular frequency E_curr. Specifically, if the medium of the integrated circuit is lossy media tri and tetra, the angular frequency under the lossy medium is ω_{curr_a}Will be ω_{curr_a}＝2πF_{min_a}Substituting the matrix equation of formula (12a) to obtain omega_{curr_a}Field solution at angular frequency E_{curr_a}(ii) a If the medium of the integrated circuit is lossless media (i) and (ii), the angular frequency under the lossless media is omega_{curr_b}Will be ω_{curr_b}＝2πF_{min_b}Substituting the matrix equation of formula (12b) to obtain omega_{curr_b}Field solution at angular frequency E_{curr_b}。

Step 423, calculating the field solution E by a relative error calculation formula corresponding to the type of dielectric loss of the integrated circuit_currWherein the relative error res is used to evaluate whether the field solution of the frequency point is accurate. Specifically, if the medium of the integrated circuit is lossy media (c) and (c), the field solution E is calculated by a relative error calculation formula of the following formula (16a)_{curr_a}Relative error res _ a of:

if the media of the integrated circuit are lossless media (i) and (ii), calculating the field solution E by a relative error calculation formula of the following formula (16b)_{curr_b}Relative error res _ b:

in the formulas (16a) and (16b), the numerator is the residual error, the denominator is the source term, the modulus of the two is the relative error of the two, and E_currFor the current angular frequency omega_currAnd (4) field solution.

Step 424, when the relative error res is less than or equal to epsilon₁When the frequency is the lowest, the description shows that the precision requirement is met when the iteration of the current round is started, and the reference frequency point f is obtained_ref＝ω_currPer 2 pi and field solution E thereof_ref＝E_currAnd ending; at the relative error res>ε₁It jumps to step 425. Wherein epsilon₁Is a preset upper limit of the error threshold. Specifically, if lossy media are three and four, reference frequency point f is obtained_ref＝ω_{curr_a}/2 π, field solution E thereof_ref＝E_{curr_a}If the media are lossless, f_ref＝ω_{curr_b}/2 π, field solution E thereof_ref＝E_{curr_b}；

Step 425, convert ω to ω_curr＝π(F_min+F_max) Substituting the matrix equation corresponding to the dielectric loss type of the integrated circuit to obtain a new field solution E_curr'. Specifically, in the case of lossy media tri and tetra, ω will be represented by_{curr_a}＝π(F_{min_a}+F_{max_a}) Substitution of formula (12a) to obtain a new field solution E under lossy media_{curr_a}'; under the condition of no damage to media (I) and (II), the value of omega is adjusted_{curr_b}＝π(F_{min_b}+F_{max_b}) Substitution of formula (12b) to obtain a new field solution E under lossless media_{curr_b}′。

Step 426, calculating the new field solution E by the relative error calculation formula_curr'a new relative error res' is obtained.

Step 427, in the new field solution E_curr'relative error res'<ε₀When making F_max＝ω_currA/2 pi and jumps to step 425; the relative error res' of the new field solution is less than or equal to epsilon₁And res' is more than or equal to epsilon₀Then, a reference frequency point f is obtained_ref＝ω_currAnd/2 pi and field solution thereof, and finishing; relative error res 'at the new field solution'>ε₁When making F_min＝ω_currAnd/2 pi and jumps to step 425. Wherein epsilon₀Is a predetermined lower limit of the error threshold, and e₀<ε₁，ε₀Can take the value of 10^-5。

Specifically, in the case of lossy media (c) and (c), the new relative error obtained in step 426 is res _ a ', in res _ a'<ε₀Then, make the new iteration frequency upper limit value F_{max_a}′＝ω_{curr_a}A value of/2 pi, and_{max_a}value of' as F_{max_a}Substituting step 425; and in res _ a'>ε₁Then, the new lower limit value F of the iteration frequency is set_{min_a}′＝ω_{curr_a}A value of/2 pi, and_{min_a}value of' as F_{min_a}Step 425 is substituted.

Similarly, in the case of lossless media (r) and (r), the new relative error obtained at step 426 is res _ b ', at res _ b'<ε₀Then, make the new iteration frequency upper limit value F_{max_b}′＝ω_{curr_b}A value of/2 pi, and_{max_b}value of' as F_{max_b}Substituting step 425; and res _ b'>ε₁Then, the new lower limit value F of the iteration frequency is set_{min_b}′＝ω_{curr_b}A value of/2 pi, and_{min_b}value of' as F_{min_b}Step 425 is substituted.

The error threshold is a quantitative index of whether the relative error meets the standard, if the error does not meet the standard, the frequency is still invalid (the field solution is inaccurate), if the error meets the standard, the frequency is reliable (the field solution is accurate), but the reliable frequency is not necessarily a reference frequency point, because the frequency may not be the lowest reliable frequency, the lowest reliable frequency can be obtained through iteration and used as the reference frequency point. This embodiment sets up ∈₀And ε₁In step 427, the two thresholds are ended as long as the error is between the two thresholds, so that the convergence rate of the bisection iteration can be accelerated, if only one threshold is set in the prior art, the error is greater than the threshold and goes left, and the error is smaller than the threshold and goes right, so that the error can be really close to the threshold by iteration for many times, and the convergence rate is slow.

It should be understood that steps 421 to 427 are sequentially executed steps, and are executed according to the step number specified when there is a jump in the step content, for example, there is an iterative process through jumps in step 427, as long as the calculated relative error of each iteration is less than epsilon₀Iterations occur and the field solution and relative error are recalculated, and so long as the relative error calculated for each iteration is greater than ε₁Iteration will also occur and the field solution and relative error will be recalculated until the relative error calculated after iteration is at ε₀And ε₁Thereby realizing the calculation of the reference frequency point through the step 421 and 427.

Step 421 and 427 can calculate the reference frequency point and the solution at the reference frequency point by using an iterative method based on the critical frequency point, so that the electromagnetic field simulation field solution at which frequency is obtained is definitely accurate.

When the dielectric is lossless dielectric, matrix K is in nanometer level in the very large scale integrated circuit₁And K₂By 10 of the element of³⁴Of the order of (1), so adding or subtracting two matrices directly at low frequencies results in a matrix omega²K₂Considered as zero. In order to solve the problem fundamentally, a generalized eigenvalue decomposition method is provided, and an eigenvector of the matrix is extracted, wherein the eigenvector is only related to the scale characteristic of the integrated circuit and is not related to the simulation frequency.

For the matrix K₁And K₂The generalized eigenvalue problem of frequency independence is shown in formula (17):

K₁x＝λK₂x (17)；

in the formula (17), λ is an eigenvalue, x is an eigenvector, and the formula (17) has N sets of eigenvalues and solution of the eigenvector in common, that is, (λ)₁,λ₂,…,λ_N) And (x)₁,x₂,…,x_N). Due to K₁Is symmetrically semi-positive, K₂Are symmetrically positive, so that the eigenvalues λ are non-negative and real, while x and K are₁、K₂Are orthogonal. In addition, the medium is a non-dispersive medium, so that the matrix K₁And K₂Is not changed with the change of the frequency, so the eigenvalues and eigenvectors obtained by the generalized eigenvalue problem of equation (17) are applicable to all frequencies.

Let Φ equal to (x)₁,x₂,…,x_N) Then obtaining formula (18) and formula (19):

Φ^TK₁Φ＝Λ (18)；

Φ^TK₂Φ＝I (19)；

where φ is the matrix formed by all eigenvectors, Λ is the diagonal matrix formed by the eigenvalues, and I is an identity matrix.

The inverse matrix of equation (20) is obtained from equations (18) and (19):

K(ω)^-1＝Φ(Λ-ω²I)^-1Φ^T (20)；

since the eigenvalues in the formula (17) can be divided into two groups, one group is related to the physical direct current mode and the non-physical direct current mode generated by the null space of the solution, the eigenvalue is theoretically zero, and the corresponding eigenvalue and eigenvector are respectively marked as Λ₀And phi₀(ii) a The other group is related to the non-zero resonance frequency of the three-dimensional structure of the integrated circuit, is a high-order eigenmode, the eigenvalue of the eigenmode is larger than zero, and the corresponding eigenvalue and eigenvector are respectively marked as lambda_highAnd phi_highAnd due to Λ₀By changing equation (20) to equation (21) when 0:

in the formula (21), phi₀The feature vector corresponding to the zero eigenvalue, called the null space vector, Φ_highFor eigenvectors corresponding to non-zero eigenvalues, Λ_highIs a diagonal matrix of non-zero eigenvalues, phi_highAnd Λ_highI.e. higher order modes associated with non-zero resonant frequencies.

As can be seen from equation (21), the frequency dependence of the solution of the matrix equation shown in equation (12b) is clearly expressed. To the right of equation (21), all other expressions are frequency independent, except for the expression explicitly containing ω. With such a frequency continuous function, a field solution of any low frequency of the integrated circuit from a high frequency to a full wave range including a direct current can be strictly obtained without the problem that the solver fails at the low frequency.

Since, as can be seen from equation (21), given an arbitrary frequency ω, the field solution is a superposition of several three-dimensional eigenmodes. For the direct current eigenmodes, i.e., eigenvectors corresponding to zero eigenvalues, their weights in the field solution are (1/ω)²) Is in direct proportion; for higher order eigenmodes, their weight in the field solution for the ith eigenmode is proportional to 1/(λ ^ k)_i-ω²) Wherein λ is_iIs the characteristic value corresponding to the ith eigenmode, and i is the serial number of the eigenmode.

At low frequencies, when the higher eigenmodes are weighted significantly less than the dc eigenmodes, i.e., [1/(λ ^ k)_i-ω²)]＜＜1/ω²The contribution of the higher order eigenmodes to the field solution is negligible. Therefore, expression (21) may become expression (22):

from this, the solution of the integrated circuit full-wave electromagnetic simulation matrix equation shown in equation (12b) below the lossless dispersion-free medium is equation (23):

as can be seen from equation (23), if the eigenvectors corresponding to all the zero eigenvalues of the generalized eigenvalue problem shown in equation (17) are solved, the field at low frequency can be accurately solved.

For a three-dimensional multi-layer VLSI, although the number of physical DC modes may be small, the null space mixes the physical DC mode and the non-physical DC mode, and the linear combination of the two modes exists in the same null space. Therefore, the physical dc mode and the non-physical dc mode cannot be distinguished from only the null space vector.

Furthermore, the size of the nullspaces cannot be reduced by discarding subsets of the nullspaces in order to increase the computation speed during the solution process, because the subsets are linearly independent of each other and each of them is essential for establishing a complete null, so that the remainder after discarding any subset of the null is incomplete, or, given an excitation vector, it can have a projection on all the null vectors, so that each null contributes to the field solution.

However, for the most advanced current very large scale integrated circuits, which have a multi-scale structure with dimensions from centimeter to nanometer, and thus the grid cells for which field division is calculated are in the order of tens of millions or even hundreds of millions, and the matrix equations shown by the equations (12a) and (12b) thus formed are in the order of tens of millions to hundreds of millions, if all the null-space vectors are solved and stored, the calculation cost is high, and a huge memory storage is required, because the null space of the matrix equations is large and grows linearly with the matrix size. Therefore, how to deal with the increased null space becomes the key to quickly and accurately solve the response of the very large scale integrated circuit at low frequency.

The method for solving the problem is proposed by utilizing the characteristic that the generalized eigenvalue problem shown in the formula (17) has a plurality of zero eigenvalues, and all the zero space vectors share the same zero eigenvalue due to the plurality of identical zero eigenvalues, although the eigenvectors of the zero space vectors are completely different.

Based on this, the present application uses the vector at the right end of equation (12a) or (12b) (the excitation vector) to reduce the dimension of the space where the field solution is located, so the right-hand source vector is always known, i.e., for a given right-hand term b (ω) to solve equation (12a) or (12b), all null-space vectors are effectively combined together to form the field solution of the matrix equation at low frequencies, as shown in equation (23). Further, the contribution of all null space vectors representing low frequencies can be represented by a vector E₀To represent that the vector covers all low frequency solutions, as shown in equation (24):

due to phi₀The eigenvector corresponding to the zero eigenvalue, whose value is independent of the frequency, is only dependent on the three-dimensional structure of the integrated circuit, so that in the case of low frequencies, the field solution E (ω) is inversely proportional to the angular frequency ω, and if the field solution at a certain low frequency can be accurately calculated, the field solution at another low frequency can be calculated according to equation (24) without strictly solving the matrix equation shown in equation (12a) or (12 b).

Therefore, the reference frequency point f is determined through the steps_refIn one aspect of the inventionThe matrix equation shown in the formula (12a) or (12b) formed at the reference frequency point can be accurately solved, and on the other hand, the reference frequency point f is adopted_refThe obtained solution vector is expressed by E₀The solution space is formed, so that the solution of the representation field at the reference frequency point needs to have an expression form as shown in the formula (24), that is, f_refThe solution of the field should be dominated by the DC eigenmode and satisfy [1/(λ)_i-ω²)]＜＜1/ω²I.e. by

The second term at the right end of this equation (21) is negligible, i.e., the contribution of the higher order eigenmodes is negligible. λ here_iIs any non-zero eigenvalue, which satisfies lambda_min≤λ_i≤λ_maxWherein λ is_minIs the minimum non-zero eigenvalue, λ_maxIs the largest non-zero eigenvalue, λ_minAnd λ_maxCorresponding to the lowest resonant frequency f of the three-dimensional multi-scale structure integrated circuit_minAnd the highest resonance frequency f_maxLowest resonance frequency f_minMaximum resonance frequency f corresponding to maximum size of VLSI_maxThen corresponding to the minimum mesh size of the VLSI after mesh generation, therefore f_maxAnd f_minThe ratio of the maximum dimension to the minimum dimension of the subdivision grid is the ratio of the maximum dimension to the minimum dimension of the VLSI, and the ratio of the maximum non-zero eigenvalue to the minimum non-zero eigenvalue is the square of the frequency or dimension ratio.

Let λ_maxAnd λ_minThe order of magnitude of the ratio is m, the critical frequency point f₀The corresponding characteristic value is recorded as lambda₀Then the selected reference frequency point f_refIt should satisfy: f. of₀<f_ref＜＜f_minSince the resonance frequency corresponds to a value whose generalized eigenvalue is not zero, λ_iSince the non-zero eigenvalue has a large value and a high resonance frequency, even the lowest resonance frequency is much larger than the reference frequency point, and the problem of the generalized eigenvalue shown in the formula (17) is not actually solved, the lowest resonance frequency f is set to be the lowest resonance frequency f_minIs an unknown quantity.

Thus, the reference frequency point f is determined_refThe principle of (a) is to ensure that the solver does not fail under the reference frequency point and ensure that the reference frequency point f_refAs far as possible below the lowest resonance frequency f_minThis ensures that the second term to the right of equation (21) is negligible, i.e., the contribution of the higher order eigenmodes is negligible.

And 500, for frequency points to be solved which are lower than the reference frequency point, acquiring field solutions under the frequency points to be solved in a mode corresponding to the medium type, wherein when the medium type belongs to a lossy medium, the field solutions of the reference frequency point are firstly partitioned, the field solutions under the frequency points to be solved are acquired according to the property of the characteristic value of the matrix of the finite element system and the partitioning result, and when the medium type belongs to a dispersive medium, the proportionality coefficient between the field solutions of the reference frequency point and the low-frequency field solutions is firstly calculated, and the field solutions under the frequency points to be solved are acquired according to the proportionality coefficient.

Four types of media, i.e., lossless and frequency-dispersion-free media, lossy and frequency-dispersion-free media, and lossy and frequency-dispersion-free media, are introduced below.

The calculation method of the frequency point field solution to be solved when the medium type is a loss-free and dispersion-free medium is as follows.

In the case of formula (i), formula (11) is substituted for formula (24) to obtain formula (25), and it can be understood that the reference frequency point is a frequency point at which the field solution can be accurately obtained by the finite element method, and also belongs to the low frequency range, that is, the field solution E (ω) at the reference frequency point satisfies formula (25):

wherein, because the medium is a lossless and non-dispersive medium, the magnetic permeability in all the units e

And dielectric constant

Are independent of frequency, so that the quantities other than ω in the right-hand term of equation (25) are constant values as the angular frequency ω varies, i.e. the field solution E (ω) is inversely proportional to the angular frequency ω in the low frequency case, i.e. at the reference frequency f_refIs accurate solution E_refWhen known, can be based on field solutions E_refAlso satisfying the properties of formula (25), gives formula (26):

where j is an imaginary unit, μ₀Is the magnetic permeability of a vacuum medium, phi₀Is the eigenvector corresponding to the zero eigenvalue, T represents the transposition of the matrix, e is the basic unit obtained by mesh division, V^eIs an integral body of the basic unit e,

for the ith interpolation basis function of the base unit e,

is the external excitation current source of the basic unit e. The lower than reference frequency f is described by equation (26)_refIs able to solve the field by E_refBy performing a reverse-derivation, in one embodiment, the solution e (f) for the frequency f to be found can be calculated by equation (27):

E(f)＝(f_refE_ref)/f (27)；

therefore, in order to calculate the electromagnetic response of the very large scale integrated circuit at the frequency point lower than the reference frequency point, the reference frequency point and the accurate field solution at the reference frequency point can be calculated first, and then the electromagnetic response of the very large scale integrated circuit at the low frequency lower than the reference frequency point can be obtained by analyzing the formula (27).

The reference frequency point f_refField solution E of_refInstead of being derived from equation (26), equation (26) is only used to demonstrate that the frequency f to be determined is calculated from equation (27)The solution is accurate, effective and feasible, and the field solution E is actually calculated_refWhat is needed is a step of simultaneously obtaining a field solution in each iteration of calculating the reference frequency points in step 420, where each iteration of one frequency point will calculate a corresponding field solution, for example, an angular frequency ω is obtained in step 425_currField solution E of_curr' if it is determined in step 427 that the angular frequency is the reference frequency point, the corresponding field solution E_curr' is the field solution E of the reference frequency point_ref。

Therefore, the scheme of the integrated circuit full-wave electromagnetic simulation under the condition of no loss and no frequency dispersion medium can be obtained: firstly, the reference frequency point of the integrated circuit to be simulated is calculated through the step 421 and the step 427, and for the frequency band higher than the reference frequency point, the solution can be directly carried out by adopting a full three-dimensional electromagnetic field numerical calculation method; for the frequency bands lower than the reference frequency point, the solution of the frequency point to be solved is obtained through the solution of the reference frequency point, so that the problem of continuity of low-frequency and high-frequency responses is solved, the simulation result is more accurate, and the problem of discontinuous response splicing curves of two solvers caused by the difference of solution results at the frequency point where the high frequency band and the low frequency band intersect when different solvers are respectively adopted for the high frequency and the low frequency is avoided.

And when the medium type is a loss-free dispersive medium, the calculation mode of the frequency point field solution to be solved is as follows.

The same as (i), in case of (ii), the same as (25) is also required, the same parts as (ii) are not described herein, and different from (ii), after obtaining (25), when the dielectric constant and the magnetic permeability of the dielectric layer forming the multi-layered lsi have dispersion characteristics, the matrix K is different from (ii)₁And K₂Is no longer frequency independent, i.e. the eigenvector Φ corresponding to the zero eigenvalue of the above equation₀The frequency is not related to the frequency any more, and at the moment, the solution of the frequency to be solved cannot be directly solved according to the solution under the reference frequency point. However, the dimension of the space in which the field solution is located can still be reduced by using the vector of the right-hand term b (ω) (excitation vector), i.e. for a given right-hand term, all low-frequency solutions can be covered by one vector, i.e. the contribution of all null-space vectors representing low frequencies can be usedA vector E₀Thus, in the case of dielectric constant and permeability of dielectric layers forming a multi-layer LSI circuit having dispersion characteristics, the low frequency solution can still use a vector E representing the contribution of all the zero-space vectors of the low frequency₀Thus, the low frequency solution has a scaling relationship with the solution of the reference frequency point, and therefore, in one embodiment, the scaling relationship is given as equation (28):

E(f)＝kE_ref (28)；

where k is a proportionality coefficient and k is a real number, E_refIs a reference frequency point f_refThe exact field solution of (c) is obtained by step 420.

Formula (29) can be obtained by substituting formula (28) for formula (12 b):

k(K₁-ω²K₂)E_ref＝b(ω) (29)；

since the matrix K is in the nanometer range of discrete tetrahedral size in the very advanced VLSI circuits₁And K₂The phase difference is of the order of 10³⁴Directly performing ω at low frequency²K₂And K₁In the subtraction of (2), the computer directly uses the matrix omega²K₂Considered to be zero, the value of k is still not solved accurately. But note that since this time E_refIs representative of a low frequency solution, thus E_refThe method is the combination of all zero space vectors representing low frequency, and the corresponding characteristic value of the zero space vector is 0, namely K is obtained₁Φ₀0, to give K₁E_ref0. Will K₁E_refSubstitution of formula (29) with 0 yields formula (30):

k(ω²K₂)E_ref＝b(ω) (30)；

both ends of the pair (30) are multiplied by non-zero vectors

To give formula (31):

the formula for calculating the scaling factor k can be obtained from the formula (31), that is, the formula (32):

as can be seen from equation (32), at low frequencies, [1/(λ) ]_i-ω²)]＜＜1/ω²When the time is short, namely the contribution of the high-order eigenmode to the field solution can be ignored, the field solution of the matrix equation under the lossless dispersive medium can pass through the reference frequency point f_refIs accurate solution E_refCalculated based on the reference frequency point f_refIs accurate solution E_refCalculating the frequency f (f) to be solved<f_ref) The following field solution is calculated by equation (33):

the above formula shows that, in order to calculate the electromagnetic response of the very large scale integrated circuit at the frequency point lower than the reference frequency point, the accurate solution at the reference frequency point can be calculated first, and then the electromagnetic response of the very large scale integrated circuit at the low frequency lower than the reference frequency point can be obtained by analyzing the formula (33).

The scheme of the full-wave electromagnetic simulation of the integrated circuit under the condition of no loss and frequency dispersion medium is similar to the scheme of the lossless and frequency dispersion medium, and the difference is that the solution of the frequency point to be solved is obtained through the field solution of the reference frequency point and the proportionality coefficient between the reference frequency point and the low frequency solution, but the problem of continuity of low-frequency and high-frequency responses can be solved, the simulation result is more accurate, and the problem of discontinuous response splicing curves of two solvers caused by difference of the solution results at the frequency point where the high frequency band and the low frequency band intersect when different solvers are adopted for the high frequency and the low frequency is avoided.

And when the medium type is a lossy non-dispersive medium, the calculation mode of the frequency point field solution to be solved is as follows.

In case of (c), before the matrix equation under the lossy dispersion-free medium is established in step 300, the following steps may be performed: and step C, dividing the material for forming the multilayer integrated circuit into a non-conductive medium and a conductive medium, numbering the grid nodes based on the divided parallel flat plate field according to the medium areas, continuously arranging the grid nodes according to the medium areas, taking the numbering and arranging results as the grid division results, and then establishing a matrix equation under the lossy non-frequency-dispersion medium based on the grid division results and the wave equation of the electric field in step 300.

Specifically, when the material forming the multilayer integrated circuit is a lossy dispersion-free medium, the medium is composed of a non-conductive medium and a conductive medium, the non-conductive medium is usually located between copper-clad layers of the multilayer integrated circuit to prevent the copper-clad layers of different layers from directly contacting, and the conductive medium is the copper-clad layer of the multilayer integrated circuit and is used for transmitting power required by different components of the integrated circuit during operation to each position of the integrated circuit and transmitting high-speed signals.

Numbering the areas where materials for forming the multilayer integrated circuit are located according to the non-conductive media and the conductive media, and establishing the mapping relation between the area numbers and the non-conductive media and the conductive media; then, based on the established mapping relation between the area numbers and the non-conductive media and the conductive media, the grid nodes of the split parallel flat plate field are numbered again according to the medium areas where the grid nodes are located, and the grid nodes are continuously arranged according to the medium areas where the grid nodes are located and the non-conductive media and the conductive media. For example, assuming that 100000 mesh nodes exist in the whole divided mesh, wherein 70000 mesh nodes are located in the non-conductive medium region, 30000 mesh nodes are located in the conductive medium region, and for a mesh node located at the interface between the conductive medium and the non-conductive medium, it is considered that the mesh node belongs to the conductive medium region, all mesh nodes located in the non-conductive medium region may be renumbered to be 1 to 70000, and all mesh nodes located in the conductive medium region may be renumbered to be 70001 to 100000.

It will be appreciated that because the mesh nodes have been numbered during the mesh generation of step B2 without regard to the media regions in which the mesh nodes are located, the numbered mesh nodes are renumbered in this step because this step is a sequential numbering of the non-conductive media and conductive media regions and therefore differs from the numbering in step B2.

In case of (c), in one embodiment, before calculating the reference frequency point and the field solution of the integrated circuit in step 400, the following steps may be performed: and step 360, dividing the finite element unknowns based on the numbering result of the grid nodes.

The finite element unknown quantity is determined according to the wave equation to be solved, if the wave equation of the electric field is to be solved, the finite element unknown quantity is the electric field of the basic unit edge or the surface element, and if the wave equation of the magnetic field is to be solved, the finite element unknown quantity is the magnetic field of the basic unit edge or the surface element. In this embodiment, the wave equation of the electric field is solved, so that the finite element unknowns refer to the electric field of the cell edge or surface element.

Step 360 may include step 361 and step 362.

And 361, dividing the finite element unknowns according to the dividing results of the medium regions according to the number of the grid nodes in the non-conductive medium regions and the numbering results of the grid nodes.

Since conductivity cannot be neglected in lossy media, the portion of the stiffness matrix that is related to conductivity (i.e., K)₃) The method becomes one of key bases for solving the field solution under the lossy medium, and the integrated circuit needs to be divided into regions according to the conductive medium and the non-conductive medium.

The division rule of the finite element unknown quantity should be the same as the numbering rule of the grid nodes, for example, if the numbering sequence of the grid nodes is that the non-conductive medium region is in front and the conductive medium region is behind, the division of the finite element unknown quantity should also be divided according to the sequence that the non-conductive medium region is in front and the conductive medium region is behind. In the dividing process, the finite element unknown quantity in the non-conductive medium region is recorded as n, the finite element unknown quantity in the conductive medium region is recorded as c, and the finite element unknown quantity at the interface of the conductive medium and the non-conductive medium is included in c.

And step 362, partitioning the stiffness matrix of the equation set according to the sequencing result of the finite element unknowns.

After the reordering and numbering, the stiffness matrix may be partitioned, and the sequence of partitioning the matrix should be the same as the partitioning sequence of the finite element unknowns, for example, the partitioning for the above finite element unknowns should also be a rule of partitioning in the sequence of the non-conductive medium region before the conductive medium region, and the left matrix K (ω) ═ K of equation (12a)₁-ω²K₂+jωK₃After blocking, this becomes:

wherein, K_nn(ω) is a submatrix associated with a region of a non-conductive medium, K_nc(ω) is a submatrix of non-conductive medium regions associated with conductive medium regions, K_cn(ω) is a submatrix of regions of conductive media associated with regions of non-conductive media, i.e., K_ncTransposition of (ω), K_cc(ω) is a submatrix associated with the region of conductive medium. The purpose of numbering in step 361 and ordering in this step is to form K_nn(ω)、K_nc(ω)、K_cn(omega) and K_cc(ω) and making the media area correspond to the index of the block matrix.

Each block submatrix in the above formula is specifically:

wherein, K_1,nnIs a matrix K₁A submatrix associated with the non-conductive media region;

K_2,nnis a matrix K₂A submatrix associated with the non-conductive media region;

K_1,ncis a matrix K₁A sub-matrix in which the medium non-conductive medium region is associated with the conductive medium region;

K_2,ncis a matrix K₂A sub-matrix in which the medium non-conductive medium region is associated with the conductive medium region;

K_1,cnis a matrix K₁A submatrix associated with the medium conductive medium region and the non-conductive medium region;

K_2,cnis a matrix K₂A submatrix associated with the medium conductive medium region and the non-conductive medium region;

K_1,ccis a matrix K₁A sub-matrix associated with the region of conductive medium;

K_2,ccis a matrix K₂A sub-matrix associated with the region of conductive medium;

K_3,ccis a matrix K₃With a sub-matrix associated with the region of conductive medium.

Due to the matrix K₃Containing only conductivity variables, thus the matrix K₃Present only in K_cc(ω) in the submatrix.

In the case of (iii), in one embodiment, the blocking the field solution of the reference frequency point in step 500 includes: and D1, partitioning the field solution of the reference frequency point based on the partitioning result of the finite element unknowns. The reference frequency point means that the field solution of the integrated circuit calculated at the reference frequency point is definitely accurate, the direct current eigenmode of the electromagnetic wave at the reference frequency point plays a main role, and the contribution of the high-order eigenmode can be ignored.

The step D1 specifically includes: partitioning the field solution of the reference frequency point according to the partitioning result of the stiffness matrix of the equation set to obtain

Wherein E is_refField solutions for reference frequency points, E_n,refFor the part of the field solution at the reference frequency lying within the non-conductive medium region, E_c,refIs the part of the field solution under the reference frequency point in the conductive medium area.

After the finite element unknowns are reordered and numbered, and the finite element sparse matrix is partitioned, when the simulated angular frequency is omega, the field solution partition at the left end of the formula (12a) is the following formula (34):

wherein E is_n(ω) is the portion of the field solution lying within the region of the non-conductive medium at an angular frequency ω, E_c(ω) is the portion of the field solution that lies within the region of the conductive medium at an angular frequency ω.

Similarly, the right-hand term of equation (12a) is also partitioned to yield: b (ω) ═ b_n(ω),b_c(ω)]^TWherein b is_n(ω) is the portion of the source term that lies within the region of the non-conductive medium, b_c(ω) is the portion of the source term that is located within the region of the conductive medium. Since it can be considered that the current source exists only in the non-conductive medium region for excitation, for example, the current density value of the excitation current region is set, the region conductivity of the excitation current is considered to be 0, and the skin effect is avoided from influencing the current distribution, the right-end term can be written as the following formula (35):

wherein b (ω) is a source term when the angular frequency is ω, b₀Is a constant vector independent of frequency, V^eIs an integral of the basic cell e.

Thus, equation (12a) is at reference frequency f_refField solution E_refIs represented by the following formula (36):

In thirdIn one embodiment, the obtaining of the to-be-solved frequency point lower field solution according to the property of the finite element system matrix eigenvalue and the blocking result in step 500 includes: step D2, splitting the field solution of the reference frequency point into a real part and an imaginary part, and setting the imaginary part in the conductive medium area to zero to obtain the split field solution of the reference frequency point:

wherein f is the frequency to be solved, f_refIs a reference frequency point.

For the frequency bands higher than the reference frequency point, a full three-dimensional electromagnetic field numerical calculation method can be directly adopted for solving; for the frequency bands lower than the reference frequency point, the frequency points to be solved are obtained by the solution of the reference frequency point, so that the problem of continuity of low-frequency and high-frequency responses is solved, the simulation result is more accurate, and the problem of discontinuous response splicing curves of two solvers caused by the difference of the solution results of the frequency points where the high-frequency band and the low-frequency band intersect when different solvers are respectively adopted for the high frequency and the low frequency is avoided.

Specifically, after the chunking in step D1 is completed, in combination with the characteristics of the right-end source term in formula (35) in step D1, that is, the right-end source term exists only in the non-conductive medium region, and the right-end source term is a pure imaginary number, for the field solution E (ω) of formula (34), the field solution E is located in the non-conductive medium region_n(omega) and a field solution E in the region of the conductive medium_c(ω) is the following formula (37):

wherein, the matrix

Φ_0,nnThe eigenvector corresponding to the zero eigenvalue in the initial generalized eigenvalue equation corresponding to the block matrix of the non-conductive medium region, that is, the eigenvector corresponding to the zero eigenvalue in the following formula (38);

Φ_high,nnthe characteristic vector corresponding to the non-zero eigenvalue in the initial generalized eigenvalue equation corresponding to the block matrix of the non-conductive medium area, namely the characteristic vector corresponding to the non-zero eigenvalue in the following formula (38);

Λ_high,nna diagonal matrix of non-zero eigenvalues in an initial generalized eigenvalue equation corresponding to the block matrix of the non-conductive medium region, that is, a diagonal matrix of non-zero eigenvalues in the following equation (38);

K_1,nnx＝λK_2,nnx (38)；

Φ_0,ccthe eigenvector corresponding to the zero eigenvalue in the initial generalized eigenvalue equation corresponding to the block matrix of the conductive medium region, that is, the eigenvector corresponding to the zero eigenvalue in the following formula (39);

wherein,

thus, the complex field solution E (ω) represented by equation (37) is divided into a real part re (E (ω)) and an imaginary part im (E (ω)) represented by equation (40):

as can be seen from equation (40), at low frequencies, the real part of the field solution is independent of frequency, while the imaginary part is partitioned into non-conductive medium regions and conductive medium regions, the field solution for the non-conductive medium regions is inversely proportional to angular frequency or frequency, and the field solution for the conductive medium regions is 0.

The reference frequency point f is obtained through the step 421 and 427_refThen, after the separation of the real part and the imaginary part of e (f) is completed through step D2, the field solutions at other frequency points lower than the reference frequency point can be calculated based on the reference frequency point according to the following formula, so as to obtain the low-frequency electromagnetic response of the lsi below the reference frequency point.

Specifically, at the reference frequency point f_refField solution E of_refAfter the splitting is completed, the real part of the signal is kept unchanged, and for the imaginary part of the signal, the signal is split into a non-conductive medium area part and a conductive medium area part, the part of the conductive medium area is set to zero, and the part of the non-conductive medium area is processed by the formula im (f)_refE_n,ref) Calculating the/f, namely finally obtaining the field solution of the frequency to be solved according to the following formula (41):

and when the medium type is the lossy dispersive medium, the calculation mode of the frequency point field solution to be solved is as follows.

The same as the third one, in the case of the fourth one, the step C is also required, that is, the material forming the multilayer integrated circuit is divided into the non-conductive medium and the conductive medium, the mesh nodes based on the divided parallel flat plate field are numbered according to the medium area, so that the mesh nodes are continuously arranged according to the medium area where the mesh nodes are located, the result of the numbering and the arranging is taken as the mesh division result, and then the step 300 is performed to establish the matrix equation under the lossy dispersive medium based on the mesh division result and the wave equation of the electric field.

Before calculating the reference frequency point and the field solution of the integrated circuit in step 400, step 360 is also required, that is, the finite element unknown quantity is divided based on the numbering result of the grid nodes, and the dividing method is also the same as step 361 and step 362, the finite element unknown quantity is divided according to the dividing result of the medium region according to the number of the grid nodes in the non-conductive medium region and the numbering result of the grid nodes, and then the stiffness matrix of the equation set is divided according to the sorting result of the finite element unknown quantity.

When the field solution of the reference frequency point is blocked, the same is performed through step D1, that is, the field solution of the reference frequency point is blocked according to the blocking result of the stiffness matrix of the equation set, so as to obtain the block result

The above-mentioned parts that are the same as those in the third embodiment are not described herein again, but different from those in the third embodiment, in the case of the fourth embodiment, the obtaining of the field solution at the frequency point to be solved includes:

wherein E is_refField solutions for reference frequency points, E_n,ref,reFor field solutions at the reference frequency point lying in the real part, E, of the non-conductive medium region_n,ref,imFor field solutions at the reference frequency point in the imaginary part, E, of the non-conducting medium region_c,ref,reThe field solution of the frequency to be solved is located in the real part of the conductive medium area under the reference frequency point, and then is calculated according to the relationship between the field solution of the frequency to be solved and the split reference frequency point field solution by the following formula: e (f) zE_refWherein z is a proportionality coefficient and z is a complex number, the proportionality coefficient z being equal to z_re+jz_imWherein z is_reIs the real part of the proportionality coefficient, z_imIs the imaginary part of the proportionality coefficient, and in the split reference frequency point field solution:

wherein the intermediate variable z of the region of the conductive medium_cAnd intermediate variable z of non-conductive medium region_nRespectively as follows:

wherein, the matrix

Φ_0,ccInitial generalized eigenvalue equations corresponding to a block matrix for a region of conductive medium

The feature vector corresponding to the mid-zero feature value,

Φ_0,nninitial generalized eigenvalue equation K corresponding to a block matrix for a region of non-conductive medium_1,nnx＝λK_2,nnThe feature vector corresponding to the zero feature value in x,

Φ_high,nninitial generalized eigenvalue equation K corresponding to a block matrix for a region of non-conductive medium_1,nnx＝λK_2,nnThe non-zero eigenvalue in x corresponds to the eigenvector,

Λ_high,nninitial generalized eigenvalue equation K corresponding to a block matrix for a region of non-conductive medium_1,nnx＝λK_2,nnThe diagonal matrix of non-zero eigenvalues in x,

b₀is a constant vector independent of frequency, omega_refIs a reference frequency point f_refAngular frequency of lower, omega_ref＝2πf_refT represents transposition;

wherein,

K_1,nnis a matrix K₁In a sub-matrix associated with a region of non-conductive medium,

K_1,ccis a matrix K₁Of the sub-matrix associated with the region of conductive medium,

K_1,cnis a matrix K₁A sub-matrix associated with the region of medium conductive medium and the region of non-conductive medium,

K_1,ncis a matrix K₁A sub-matrix of medium non-conductive media areas associated with conductive media areas,

K_2,ncis a matrix K₂A sub-matrix of medium non-conductive media areas associated with conductive media areas,

K_2,nnis a matrix K₂In a sub-matrix associated with a region of non-conductive medium,

K_3,ccis a matrix K₃A sub-matrix associated with the region of conductive medium;

specifically, if the medium employed by the integrated circuit is lossy and has a frequency dispersion characteristic, the finite element system matrix is no longer independent of frequency, but is a function of frequency, and thus, the matrix a in equation (37) is also a function of frequency, and thus, the field solution and the angular frequency ω in equation (40) are not simply an explicit relationship because the matrix a therein is also a function of the angular frequency ω. At the reference frequency point, the field solution can be written as follows:

wherein E is_n,ref,reIs E_n,refReal part of (E)_n,ref,imIs E_n,refImaginary part of, E_c,ref,reIs E_c,refThe real part of (a).

Substituting re (E (ω)) and im (E (ω)) into formula (42) to obtain formula (43):

wherein, ω is_refIs a reference frequency point f_refAngular frequency of lower, omega_ref＝2πf_refAnd providing an intermediate variable z of the region of the conductive medium_cAnd intermediate variable z of non-conductive medium region_nComprises the following steps:

then the following results are obtained:

this gives:

and in the case of (iv), since the medium used by the integrated circuit is lossy and has a dispersion characteristic, the field solution at the reference frequency point represented by equation (43) is not simply an explicit relationship with the reference frequency point, and therefore the field solution at the frequency point to be solved lower than the reference frequency point cannot be directly inferred according to equation (43), however, since the dc component still plays a major role at the reference frequency point and the frequencies below the reference frequency point, there is a scaling relationship between the field solution at the frequency point to be solved lower than the reference frequency point and the field solution at the reference frequency point, and the scaling relationship can be expressed by the following equation:

E(f)＝zE_ref

wherein z is a proportionality coefficient and z is a complex number, z ═ z_re+jz_im. Converting E (f) to zE_refFormula (44) is obtained by substituting formula (12 a):

z(K₁-ω²K₂+jωK₃)E_ref＝-b(ω) (44)；

but due to the solution E of the reference frequency point_refAlso represents a low frequency solution, while knowing E_refIs the combination of all the null space vectors representing low frequencies, and the corresponding eigenvalue of the null space vector is 0, so K₁Φ₀0, thus K₁E_refBy this, formula (44) can be simplified to formula (45):

z(-ω²K₂+jωK₃)E_ref＝-b(ω) (45)；

multiplying both ends of equation (45) simultaneously by

Formula (46) can be obtained:

wherein,

is E_refThe transposed matrix of (2).

However, even if the matrix K is eliminated₁Directly calculating- ω in the formula (45)²K₂+jωK₃An accurate solution for z cannot be obtained because the matrix K is now₂And K₃The elements of (a) still have different orders of magnitude, resulting in the computer directly converting the matrix omega in the matrix operation²K₂Considered as zero. For this purpose, the left matrix K in equation (45) is further mapped₂And K₃The partitioning is performed, and the left field solution in equation (45) is split and the real and imaginary forms, i.e., equation (43), are obtained.

Based on the finite element stiffness matrix block result of step 362, it can be known that:

wherein, K_2,nnIs a matrix K₂A submatrix associated with the non-conductive media region;

K_2,cnis a matrix K₂Middle conductive medium region and non-conductive medium regionA domain-associated sub-matrix;

In one embodiment, z_reAnd z_imIs given by the following formula:

wherein,

matrix coefficients that are the real part of the scaling coefficient z,

matrix coefficients being the imaginary part of the scaling coefficient z,

is E_n,ref,imThe transpose matrix of (a) is,

is E_n,ref,reThe transposed matrix of (1), the

And

obtained by the following formula:

where ω is the angular frequency of the electromagnetic wave, E_n,ref,reIs E_n,refReal part of (E)_n,ref,imIs E_n,refImaginary part of, E_c,ref,reIs E_c,refThe real part of (a); k_2,nnIs a matrix K₂Submatrix, K, associated with areas of non-conducting medium_3,ccIs a matrix K₃With a sub-matrix associated with the region of conductive medium.

Specifically, the splitting result of the field solution of the reference frequency point is substituted into formula (46), so as to obtain formula (47):

using the property of the combination of the solution to zero space vectors at low frequencies, we can obtain:

K_2,nnAΦ_0,ccz_c-K_2,ncΦ_0,ccz_c＝0；

K_2,nnA-K_2,nc＝0；

A^TK_2,nn-K_2,cn＝0；

thus, equation (47) is simplified to the following equation:

the simplified result is substituted into formula (46) to obtain formula (48):

converting equation (48) to a matrix form of real and imaginary parts, resulting in equation (49):

wherein,

matrix coefficients that are the real part of the scaling coefficient z,

matrix coefficients being the imaginary part of the scaling coefficient z, b₀Is a constant vector independent of frequency, and:

where ω is the electromagnetic angular frequency.

It can be seen that the matrix K shown in the formula (27) no longer appears in the formula (49)₂And K₃Direct addition expression, thereby avoiding the computer directly combining the matrix K in matrix operation₂The problem is treated as zero.

And step 600, obtaining the electromagnetic response of the full frequency band based on the field solutions of all the frequency points to be obtained.

After the field solutions of the frequency points in the low-frequency band and the high-frequency band are obtained, post-processing operation can be performed based on the discrete field quantity and the calculation information required by the user, and the method specifically includes at least one of the following implementation items:

1. calculating to obtain at least one distribution item of potential distribution, current distribution, power consumption distribution, loss distribution and heat distribution of different layers of the integrated circuit power supply system under the approximate direct current or low frequency based on the calculated discrete field quantity, and drawing a chart of the obtained distribution item;

2. further calculating user-provided port parameters based on the calculated discrete field quantities, including full-band frequency response characteristics of at least one of multi-port S parameters, multi-port impedance matrices, and equivalent circuit model parameters of the integrated circuit;

3. electromagnetic radiation of the integrated circuit and/or electromagnetic interference at a location of the integrated circuit at different frequencies is further calculated based on the calculated discrete field quantities.

4. The behavioral characteristics of the components of the integrated circuit, which may be voltage-current characteristics or other characteristics/characteristics, are further calculated based on the calculated discrete field quantities, thereby further extracting the IBIS model of the integrated circuit.

After the calculation of the item to be calculated, an integrated circuit full-wave electromagnetic simulation for the respective media type is thus completed.

The significance of electromagnetic simulation of very large scale integrated circuits is that the integrated circuits are the physical realization of very complex function very large scale circuits, and from schematic diagrams of very large scale circuits to complex integrated circuit layouts, the realization ways are infinite. When designing the schematic diagram of the super-large scale circuit, only the connection relation, the logic relation and the working time sequence of the schematic diagram are considered, and the influence of the electromagnetic interference between lines on signal transmission is not considered.

However, theories and practices prove that when high-frequency alternating current is transmitted through a physical transmission line, electromagnetic waves exist around the transmission line, and the electromagnetic waves affect the transmission line around the transmission line in an electromagnetic induction mode to form electromagnetic interference. If the effect of this electromagnetic interference is much less than the signal of the transmission line itself, then this electromagnetic interference is not sufficient to affect the operation of the integrated circuit; however, if the influence of the electromagnetic interference is equal to or even greater than the signal transmitted by the transmission line, the interference is superimposed on the signal of the transmission line and transmitted as a transmission signal, thereby changing the originally transmitted signal and destroying the operation of the integrated circuit.

In addition, when designing a schematic diagram of a very large scale circuit, resistance of a circuit wire itself is not considered, and generally, all connected wires are directly considered to be ideal conductors, and voltage drop through the conductors is 0. The voltage drop will also have a fatal influence on the signal transmission of the integrated circuit, because the signal transmission of the integrated circuit actually transmits 0 and 1 signals, the 0 and 1 signals are realized by the jump of high and low levels, the integrated circuit component identifies that the high and low levels are judged by the threshold value of the level, when the level input to the component is higher than a certain reference level (threshold value), the transmitted level is judged to be high, and when the level is lower than the reference level, the transmitted level is judged to be low. Since the high and low levels of transmission are applied to the reference voltage of the transmission line, and the voltage drop on the trace is an important factor causing the instability of the reference voltage, the instability of the reference voltage directly causes errors of the transmission result.

Therefore, calculating the electromagnetic field and the voltage and current distribution at all positions on the integrated circuit layout is an important means for analyzing the voltage drop of the integrated circuit layout through simulation calculation of the electromagnetic field. In addition, it is also an important task to observe the distribution of current on the integrated circuit layout, because when the current distribution on the integrated circuit layout is too high, the slender routing is instantly burnt, so that the integrated circuit is permanently burnt out, or a certain area of the integrated circuit is heated because of too large current, and the electromagnetic transmission of materials is changed because of temperature rise of the integrated circuit board under long-time work, so that the integrated circuit cannot normally work, or even the integrated circuit is deformed and even burnt out due to long-time accumulated heating.

In the aspect of full-wave electromagnetic simulation, with the development of the design of a large-scale integrated circuit to higher frequency and the increase of circuit complexity, simulation errors caused by high-order propagation modes are generally ignored for the simulation of a high-frequency electromagnetic field. In addition, in the traditional mode, when an equivalent circuit method is adopted for analysis, equivalent elements such as capacitance and inductance do not consider the change of the elements along with the frequency, so that errors are caused. For example, for a microstrip line or an interconnection line structure, due to the characteristics of complex intersection, step, bend, open circuit, slot, etc. of the microstrip line or the stripline, signal delay, distortion, reflection, etc. effects, etc. are generated on the microstrip line or the interconnection line, and crosstalk between adjacent lines occurs, when electromagnetic waves are transmitted in multiple modes. For this reason, full-wave electromagnetic simulation techniques are typically required to analyze the integrated circuit to obtain accurate S-parameters in discontinuous mode, since full-wave analysis takes into account all possible field components and boundary conditions.

An embodiment of the integrated circuit full-wave electromagnetic simulation system based on different dielectric characteristics disclosed in the present application is described in detail below with reference to fig. 2. The present embodiment is a system for implementing the aforementioned embodiment of the integrated circuit full-wave electromagnetic simulation method.

As shown in fig. 2, the system disclosed in this embodiment mainly includes: the device comprises a medium information acquisition module, a medium type judgment module, an equation set construction module, a reference frequency point acquisition module and a to-be-solved field solution calculation module.

The medium information acquisition module is used for acquiring layer medium material information of the integrated circuit and determining medium characteristics of the integrated circuit, which relate to electric conductivity, magnetic permeability and dielectric constant, based on the layer medium material information.

The medium type judging module is used for judging the medium type of the integrated circuit about loss and dispersion characteristics based on the medium characteristics.

And the equation set building module is used for building a matrix equation corresponding to the type of the medium based on a mesh division result obtained by carrying out mesh division on the parallel flat plate field of the integrated circuit.

And the reference frequency point acquisition module is used for calculating the reference frequency point and the field solution of the integrated circuit by an iteration method based on the matrix equation, wherein the integrated circuit field solution corresponding to the reference frequency point is an accurate solution.

The to-be-solved field solution calculation module is used for acquiring field solutions under the to-be-solved frequency points for the to-be-solved frequency points lower than the reference frequency point in a mode corresponding to the medium types, wherein when the medium types belong to lossy media, the field solutions of the reference frequency points are firstly partitioned, the field solutions under the to-be-solved frequency points are acquired according to the properties of matrix characteristic values of a finite element system and partitioning results, and when the medium types belong to dispersive media, the proportionality coefficients between the field solutions of the reference frequency points and the low-frequency field solutions are firstly calculated, and the field solutions under the to-be-solved frequency points are acquired according to the proportionality coefficients.

The electromagnetic response acquisition module is used for acquiring the electromagnetic response of the full frequency band based on the field solutions of all the frequency points to be solved.

In one embodiment, the reference frequency point obtaining module calculates the reference frequency point of the integrated circuit by:

step A2, converting the current angular frequency omega_curr＝2πF_minSubstituting into a matrix equation corresponding to the integrated circuit dielectric loss type, and solving the matrix equation to obtain omega_currField solution at angular frequency E_curr(ii) a Wherein, the matrix equation corresponding to the loss type of the medium is: (K)₁-ω²K₂) E ═ b (ω), and when the dielectric loss type is a lossy medium, the corresponding matrix equation is: (K)₁-ω²K₂+jωK₃) E ═ b (ω), where ω is the electromagnetic angular frequency, j is the imaginary unit, E is the electric field, b (ω) is the external excitation source of the whole finite element system, K is₁Is the stiffness matrix of the entire finite element system, K₂Is a dielectric constant-dependent mass matrix, K, of the entire finite element system₃Is a mass matrix related to the electrical conductivity of the whole finite element system;

the formula for the relative error res _ a for lossy media is:

wherein, ω is_{curr_a}For lossy mediaAngular frequency of lower, omega_{curr_b}At angular frequency in a lossless medium, E_{curr_a}For field solutions in lossy media, E_{curr_b}Field decomposition under lossless medium;

In one embodiment, when the medium type is a lossless non-dispersive medium, the to-be-solved field solution calculation module performs field solution E according to a reference frequency point_refCalculating a field solution e (f) for the frequency f to be solved by:

wherein f is<f_ref；

When the medium type is a lossless dispersive medium, the to-be-solved field solution calculation module obtains a solution E (f) at a to-be-solved frequency point by the following formula: e (f) ═ kE_refWhere k is a proportionality coefficient and k is a real number, E_refIs a reference frequency pointf_refThe following accurate field solutions, the proportionality coefficient k being:

wherein,

is a non-zero vector.

In one embodiment, when the media type is lossy non-dispersive media or lossy dispersive media:

In one embodiment, when the media type is lossy, frequency-dispersion free media: the to-be-solved field solution calculation module acquires a field solution under a to-be-solved frequency point through the following steps:

z is_reAnd z_imIs given by the following formula:

wherein,

matrix coefficients that are the real part of the scaling coefficient z,

matrix coefficients being the imaginary part of the scaling coefficient z,

is E_n,ref,imThe transpose matrix of (a) is,

the above-mentioned

And

obtained by the following formula:

wherein, K_2,nnIs a momentMatrix K₂Submatrix, K, associated with areas of non-conducting medium_3,ccIs a matrix K₃With a sub-matrix associated with the region of conductive medium.

The above description is only for the specific embodiments of the present application, but the scope of the present application is not limited thereto, and any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope of the present application should be covered within the scope of the present application. Therefore, the protection scope of the present application shall be subject to the protection scope of the claims.

Claims

1. An integrated circuit full-wave electromagnetic simulation method based on different dielectric characteristics is characterized by comprising the following steps:

2. The method for electromagnetic simulation of the full wave of the integrated circuit based on different dielectric characteristics as claimed in claim 1, wherein the calculating the reference frequency point of the integrated circuit by an iterative method based on the matrix equation comprises:

the formula for the relative error res _ a for lossy media is:

step A4, when the relative error res is less than or equal to epsilon₁Then, a reference frequency point f is obtained_ref＝ω_currPer 2 pi and field solution E thereof_ref＝E_currAnd ending; at the relative error res > epsilon₁Jumping to step a 5;

step A7, the relative error at the new field solution is less than ε₀When making F_max＝ω_currA/2 pi and jumping to step A5; the relative error at the new field solution is less than or equal to epsilon₁And is greater than or equal to epsilon₀Then, a reference frequency point f is obtained_ref＝ω_currAnd/2 pi and field solution thereof, and finishing; the relative error at the new field solution is greater than epsilon₁When making F_min＝ω_currA/2 pi and jumping to step A5; wherein epsilon₀Is a preset lower error threshold value epsilon₁Is a preset upper limit of the error threshold, epsilon₀＜ε₁。

3. The method for full-wave electromagnetic simulation of an integrated circuit based on different dielectric characteristics as claimed in claim 2, wherein when the dielectric type is lossless and dispersion-free dielectric, the electromagnetic simulation is performed according to a field solution E at a reference frequency point_refCalculating a field solution e (f) for the frequency f to be solved by:

wherein f is less than f_ref；

When the medium type is a lossless dispersive medium, obtaining a solution E (f) at a frequency point to be solved by the following formula: e (f) ═ kE_refWhere k is a proportionality coefficient and k is a real number, E_refIs a reference frequency pointf_refThe following accurate field solutions, the proportionality coefficient k being:

wherein,

is a non-zero vector.

4. The method for full-wave electromagnetic simulation of an integrated circuit based on different dielectric properties of claim 3, wherein when the type of medium is lossy non-dispersive medium or lossy dispersive medium:

Wherein n is a finite element unknown quantity in the non-conductive medium region, c is a finite element unknown quantity in the conductive medium region, and the finite element unknown quantity is an electric field of the basic unit edge or surface element; k (omega) is the left-end matrix of the finite element system matrix, K_nn(ω) is a submatrix associated with a region of a non-conductive medium, K_nc(ω) is a submatrix of non-conductive medium regions associated with conductive medium regions, K_cn(omega) is a region of a conductive medium and a non-conductiveSubmatrix associated with an area of an electric medium, K_cc(ω) is a submatrix associated with a region of conductive medium; and,

Wherein E is_n，refFor the part of the field solution at the reference frequency lying within the non-conductive medium region, E_c，refIs the part of the field solution under the reference frequency point in the conductive medium area.

5. The method for full-wave electromagnetic simulation of an integrated circuit based on different dielectric properties of claim 4, wherein when the type of media is lossy, frequency-dispersion free media: the acquisition of the field solution under the frequency point to be solved comprises the following steps:

wherein E is_n，ref，reFor field solutions at the reference frequency point lying in the real part, E, of the non-conductive medium region_n，ref，imFor field solutions at the reference frequency point in the imaginary part, E, of the non-conducting medium region_c，ref，reThe field solution of the frequency to be solved is located in the real part of the conductive medium area under the reference frequency point, and then is calculated according to the relationship between the field solution of the frequency to be solved and the split reference frequency point field solution by the following formula: e (f) zE_refWherein z is a proportionality coefficient and z is a complex number, the proportionality coefficient z being equal to z_re+jz_imWherein z is_reIs the real part of the proportionality coefficient, z_imIs the imaginary part of the scaling factor;

z is_reAnd z_imIs given by the following formula:

wherein,

matrix coefficients that are the real part of the scaling coefficient z,

matrix coefficients being the imaginary part of the scaling coefficient z,

is E_n，ref，imThe transpose matrix of (a) is,

is E_n，ref，reTransposed matrix of b₀Is a constant vector independent of frequency;

the above-mentioned

And

obtained by the following formula:

wherein, K_2，nnIs a matrix K₂Submatrix, K, associated with areas of non-conducting medium_3，ccIs a matrix K₃With a sub-matrix associated with the region of conductive medium.

6. An integrated circuit full-wave electromagnetic simulation system based on different dielectric characteristics, comprising:

7. The integrated circuit full-wave electromagnetic simulation system based on different medium characteristics according to claim 6, wherein the reference frequency point obtaining module calculates the reference frequency point of the integrated circuit by:

step A3, calculating the relative error res of the field solution Ecurr through a relative error calculation formula corresponding to the type of the integrated circuit dielectric loss; wherein, the calculation formula of the relative error res _ b under the lossless medium is as follows:

the formula for the relative error res _ a for lossy media is:

8. The integrated circuit full-wave electromagnetic simulation system based on different medium characteristics of claim 7, wherein when the medium type is a lossless non-dispersive medium, the to-be-solved field solution calculation module performs field solution E according to a reference frequency point_refCalculating a field solution e (f) for the frequency f to be solved by:

wherein f <f_ref；

wherein,

is a non-zero vector.

9. The integrated circuit full wave electromagnetic simulation system based on different dielectric properties of claim 8, wherein when the type of media is lossy non-dispersive media or lossy dispersive media:

Where n is a finite element unknown in the non-conductive medium region and c is a finite element unknown in the conductive medium regionA finite element unknown quantity, wherein the finite element unknown quantity is an electric field of a basic unit edge or surface element; k (omega) is the left-end matrix of the finite element system matrix, K_nn(ω) is a submatrix associated with a region of a non-conductive medium, K_nc(ω) is a submatrix of non-conductive medium regions associated with conductive medium regions, K_cn(ω) is a submatrix of regions of conductive medium associated with regions of non-conductive medium, K_cc(ω) is a submatrix associated with a region of conductive medium; and,

10. The integrated circuit full-wave electromagnetic simulation system based on different dielectric characteristics of claim 9, wherein when the type of media is lossy, frequency-dispersion free media: the to-be-solved field solution calculation module acquires a field solution under a to-be-solved frequency point through the following steps: