ACCURATE AND EFFICIENT NON LINEAR MODEL ORDER REDUCTION FOR ELECTRO-THERMAL ANALYSIS

Abstract

A method of performing an electro-thermo simulation includes defining a non-linear heat diffusion problem for at least a portion of a semiconductor device to be modeled, performing a finite volume discretization of the non-linear heat diffusion problem, reformulating a non-linear term of the discretized non-linear heat diffusion problem to decrease dimensions thereof, performing a hyper reduction of the reformulated non-linear term, and recovering the non-linear heat diffusion problem for the portion of the semiconductor device, and manufacturing the modeled semiconductor device.

Claims

1. A method, comprising: defining a non-linear heat diffusion problem for at least a portion of a semiconductor device to be modeled; performing a finite volume discretization of the non-linear heat diffusion problem; reformulating a non-linear term of the discretized non-linear heat diffusion problem to decrease dimensions thereof; performing a hyper reduction of the reformulated non-linear term, and recovering the non-linear heat diffusion problem for the portion of the semiconductor device; and causing manufacture of the modeled semiconductor device.

2. The method of claim 1, wherein performing the finite volume discretization includes defining interface nodes that represent interfaces between non-homogenous materials in the semiconductor device to be modeled.

3. The method of claim 2, wherein performing the finite volume discretization further comprises: introducing an auxiliary variable to the non-linear heat diffusion problem to redefine the non-linear heat diffusion problem to contain only quadratic non-linearities; splitting the portion of the semiconductor device to be modeled into a plurality of tetrahedrons; collocating inner thermal nodes in centroids of each of the plurality of tetrahedrons; collocating the interface nodes in centers of each contact surface between two adjacent tetrahedrons; and deriving an interconnection structure from an electrical problem associated with the non-linear heat diffusion problem.

4. The method of claim 3, wherein deriving the interconnection structure comprises defining thermal resistance values between each inner thermal node and its adjacent interface node or nodes.

5. The method of claim 4, wherein performing the finite volume discretization mathematically yields: $C\frac{\partial }{\partial t}\theta \left(t\right)+K\theta \left(t\right)+\Delta K\left(\theta \left(t\right).Math.\lambda \left(t\right)\right)=\mathrm{GP}\left(t\right)T_{0}N\frac{\partial \lambda }{\partial t}\left(t\right)+\Delta N\left(\theta \left(t\right).Math.\frac{\partial \lambda }{\partial t}\right)=M\frac{\partial \theta }{\partial t}\left(t\right)+\Delta M\left(\frac{\partial \theta }{\partial t}\left(t\right).Math.\lambda \left(t\right)\right)$ wherein C is a matrix that describes thermal capacity of the portion of the semiconductor device; wherein K is a matrix that describes thermal conductivity between each of the plurality of tetrahedrons; wherein G is a matrix that describes the power dissipation in each tetrahedron, each described powerful dissipation representing a fraction of a total power dissipation within the portion of the semiconductor device; wherein θ(t) represents unknown temperature rises in the portion of the semiconductor device; and wherein λ(t) is the auxiliary variable.

6. The method of claim 5, wherein C is a diagonal matrix with each entry c.sub.ii being nonzero provided that the index i refers to an inner thermal node, with i ranging from 1 to m.sub.e; wherein K is a matrix having a size of m.sub.e by m.sub.e, with k.sub.ij being equal to −r.sub.i if i is not equal to j, with r.sub.i being the resistance between a given inner thermal node and an interface node, and k.sub.ij being a sum in absolute value of elements of the row ij if i=j; wherein G is a diagonal matrix with each entry g.sub.ii being nonzero provided that the index i refers to an inner thermal node, with i ranging from 1 to m.sub.e; wherein N is an identity matrix having a size of m by m; and wherein M is a diagonal matrix having a size of m by m where m.sub.ij=σ(N.sub.i), with σ being a material dependent function.

7. The method of claim 6, wherein the non-linear term has a size of m.sub.e by a product of m and m.sub.e) wherein the non-linear term is reformulated by determining a squared matrix having a size of m.sub.e by m.sub.e.

8. The method of claim 7, wherein the non-linear term of the discretized non-linear heat diffusion problem is given by ΔK(θ(t).Math.λ(t)), with ΔK being non-linear; wherein the non-linear term is reformulated as: ΔK(θ(t).Math.λ(t))=A.sub.λ(t)θ(t), with A.sub.λ(t)θ(t) being the squared matrix.

9. The method of claim 8, wherein A.sub.λ(t)θ(t) is computed by replacing each diagonal value kii of the matrix K with a non-linear resistance value of λ(N.sub.it)k.sub.ii; and wherein each term λ(N.sub.i,t) is computed for the inner thermal nodes as: $\lambda \left(N_i,t\right)={\left(1+\frac{\theta \left(N_i,t\right)}{T_0}\right)}^{\sigma \left(N_i\right)}-1.$

10. The method of claim 9, wherein the hyper reduction of the reformulated non-linear term is performed by computing projection matrices able to derive a compact form for the discretized non-linear heat diffusion problem.

11. The method of claim 10, wherein the projection matrices are computed using non-linear model order reduction based upon moment matching and Krylov subspaces.

12. The method of claim 11, wherein the compact form for the discretized non-linear heat diffusion problem is: $\hat{C}\frac{\partial }{\partial t}\hat{\theta }\left(t\right)+\hat{K}\hat{\theta }\left(t\right)+\Delta \widehat{K}\left(\hat{\theta }\left(t\right).Math.\hat{\lambda }\left(t\right)\right)=\hat{G}P\left(t\right)T_0\hat{N}\frac{\partial \hat{\lambda }}{\partial t}\left(t\right)+\Delta \hat{N}\left(\hat{\theta }\left(t\right).Math.\frac{\partial \hat{\lambda }}{\partial t}\right)=\hat{M}\frac{\partial \hat{\theta }}{\partial t}\left(t\right)+\Delta \hat{M}\left(\frac{\partial \hat{\theta }}{\partial t}.Math.\hat{\lambda }\left(t\right)\right).$

13. The method of claim 12, wherein the projection matrices are defined as a matrix V having a size of {circumflex over (m)}×m.sub.e and a matrix W having a size of {circumflex over (m)}×m, where m.sub.e>m»{circumflex over (m)}.

14. The method of claim 13, wherein the hyper reduction of the reformulated non-linear term is further performed by defining reduced matrices as follows:
Ĉ=V.sup.TCV {circumflex over (m)}×{circumflex over (m)}
{circumflex over (K)}=V.sup.TKT {circumflex over (m)}×{circumflex over (m)}
Δ{circumflex over (K)}=V.sup.TΔK(V.Math.W) {circumflex over (m)}×{circumflex over (m)}.sup.2
Ĝ=V.sup.TG {circumflex over (m)}×{circumflex over (m)}
{circumflex over (N)}=W.sup.TNW {circumflex over (m)}×{circumflex over (m)}
{circumflex over (M)}=W.sup.TMV {circumflex over (m)}×{circumflex over (m)}
Δ{circumflex over (N)}=W.sup.TΔN(V.Math.W) {circumflex over (m)}×{circumflex over (m)}.sup.2
Δ{circumflex over (M)}=W.sup.TΔM(V.Math.W) {circumflex over (m)}×{circumflex over (m)}.sup.2

15. The method of claim 14, wherein ΔK is calculated as: ΔK=Σ.sub.i=1.sup.mA.sub.e.sub.i.Math.e.sub.i.sup.T; wherein ΔM is calculated as: ΔM=Σ.sub.i=1.sup.mM.Math.e.sub.i.sup.T; and wherein ΔN is calculated as: ΔN=Σ.sub.i=1.sup.mN.Math.e.sub.i.sup.T, where e.sub.i=(0, . . . , 0, 1, 0, . . . , 0).sup.T is the ith vector of a basis of an m-dimensional Euclidean space.

16. The method of claim 15, wherein the reformulated non-linear matrix ΔK has a size of {circumflex over (m)}×{circumflex over (m)}.sup.2; and wherein the hyper reduction of the reformulated non-linear term is further performed by: applying energy-conserving sampling and weighting to the reformulated non-linear matrix ΔK, beginning with a collection of a snapshots of the solution to the non-linear heat diffusion problem in a set of complex frequencies in the Laplace domain, resulting in a matrix U having dimensions of m.sub.e×s and a matrix ∧ having dimensions of m×s.

17. The method of claim 16, wherein the hyper reduction of the reformulated non-linear term is further performed by construction of a matrix H and a vector b, as follows:
H.sub.ie=V.sup.T*∧(e,i)*A.sub.pe*U(:,i),
i=1, . . . , s
e=1, . . . , m
b=H*(1, . . . , 1).sup.T; and wherein resolution of a resulting sparse minimization problem is performed as:
ξ*=argmin.sub.ξ∥Hξ−b∥.sub.F s.t.ξ≥0 and ∥ξ∥.sub.0<w wherein a solution of the problem ξ* represents a sparse vector of weights, and a set E can defined as the indices of the non-zero entries of ξ*.

18. The method of claim 17, wherein the hyper reduction of the reformulated non-linear term is further performed by evaluating it on volumed indexed by E and weighted by ξ*, as follows:
Δ{circumflex over (K)}({circumflex over (θ)}(t).Math.{circumflex over (λ)}(t))≈(Σ.sub.e=1.sup.wξ*(e)*W(e,:)*{circumflex over (λ)}(t)*V.sup.TA.sub.pe*V)*{circumflex over (θ)}(t).

19. The method of claim 18, wherein the hyper reduction of the reformulated non-linear term is further performed by performed by a singular value decomposition on the terms V.sup.T* A.sub.pe*V for e=1, . . . , w.

20. The method of claim 19, wherein recovering the non-linear heat diffusion problem for the portion of the semiconductor device is performed by solving the hyper reduction of the reformulated non-linear term to produce a reduced solution of {circumflex over (θ)}(t) and {circumflex over (λ)}(t).

21. The method of claim 20, wherein recovering the non-linear heat diffusion problem for the portion of the semiconductor device is further performed by using the projection matrices V and W as:
θ(t)=V{circumflex over (θ)}(t) λ(t)=W{circumflex over (λ)}(t).

22. A method, comprising: defining a non-linear heat diffusion problem for at least a portion of a semiconductor device to be modeled; performing a finite volume discretization of the non-linear heat diffusion problem by: defining interface nodes that represent interfaces between non-homogenous materials in the semiconductor device to be modeled introducing an auxiliary variable to the non-linear heat diffusion problem to redefine the non-linear heat diffusion problem to contain only quadratic non-linearities; splitting the portion of the semiconductor device to be modeled into a plurality of tetrahedrons; collocating inner thermal nodes in centroids of each of the plurality of tetrahedrons; collocating the interface nodes in centers of each contact surface between two adjacent tetrahedrons; and deriving an interconnection structure from an electrical problem associated with the non-linear heat diffusion problem by defining thermal resistance values between each inner thermal node and its adjacent interface node or nodes. reformulating a non-linear term of the discretized non-linear heat diffusion problem to decrease dimensions thereof; performing a hyper reduction of the reformulated non-linear term by computing projection matrices able to derive a compact form for the discretized non-linear heat diffusion problem; and recovering the non-linear heat diffusion problem for the portion of the semiconductor device from the compact form for the discretized non-linear heat diffusion problem.

23. The method of claim 22, wherein performing the finite volume discretization mathematically yields: $C\frac{\partial }{\partial t}\theta \left(t\right)+K\theta \left(t\right)+\Delta K\left(\theta \left(t\right).Math.\lambda \left(t\right)\right)=\mathrm{GP}\left(t\right)T_{0}N\frac{\partial \lambda }{\partial t}\left(t\right)+\Delta N\left(\theta \left(t\right).Math.\frac{\partial \lambda }{\partial t}\right)=M\frac{\partial \theta }{\partial t}\left(t\right)+\Delta M\left(\frac{\partial \theta }{\partial t}\left(t\right).Math.\lambda \left(t\right)\right),$ wherein C is a matrix that describes thermal capacity of the portion of the semiconductor device; wherein K is a matrix that describes thermal conductivity between each of the plurality of tetrahedrons; wherein G is a matrix that describes the power dissipation in each tetrahedron, each described powerful dissipation representing a fraction of a total power dissipation within the portion of the semiconductor device; wherein θ(t) represents unknown temperature rises in the portion of the semiconductor device; and wherein λ(t) is the auxiliary variable.

24. The method of claim 23, wherein C is a diagonal matrix with each entry c.sub.ii being nonzero provided that the index i refers to an inner thermal node, with i ranging from 1 to m.sub.e; wherein K is a matrix having a size of m.sub.e by m.sub.e, with k.sub.ij being equal to −r.sub.i if i is not equal to j, with r.sub.i being the resistance between a given inner thermal node and an interface node, and k.sub.ij being a sum in absolute value of elements of the row ij if i=j; wherein G is a diagonal matrix with each entry g.sub.ii being nonzero provided that the index i refers to an inner thermal node, with i ranging from 1 to m.sub.e; wherein N is an identity matrix having a size of m by m; and wherein M is a diagonal matrix having a size of m by m where m.sub.ij=σ(N.sub.i), with σ being a material dependent function.

25. The method of claim 24, wherein the non-linear term has a size of m.sub.e by a product of m and m.sub.e) wherein the non-linear term is reformulated by determining a squared matrix having a size of m.sub.e by m.sub.e.

26. The method of claim 25, wherein the non-linear term of the discretized non-linear heat diffusion problem is given by ΔK(θ(t).Math.λ(t)), with ΔK being non-linear; wherein the non-linear term is reformulated as: ΔK(θ(t).Math.λ(t))=A.sub.λt)θ(t), with A.sub.λ(t)θ(t) being the squared matrix.

27. The method of claim 26, wherein A.sub.λ(t)θ(t) is computed by replacing each diagonal value kii of the matrix K with a non-linear resistance value of λ(N.sub.i,t)k.sub.ii; and wherein each term λ(N.sub.i,t) is computed for the inner thermal nodes as: $\lambda \left(N_i,t\right)={\left(1+\frac{\theta \left(N_i,t\right)}{T_0}\right)}^{\sigma \left(N_i\right)}-1.$

28. The method of claim 27, wherein the projection matrices are computed using non-linear model order reduction based upon moment matching and Krylov subspaces.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] FIG. 1 is a block diagram of a computing device on which the electro-thermal simulation technique disclosed herein can be performed.

[0016] FIG. 2 is a flowchart illustrating the steps of the electro-thermal simulation technique disclosed herein.

[0017] FIG. 3 is a diagram showing sample mesh elements i and j and the contact interface therebetween, and how the electro-thermal simulation technique disclosed herein collocates a thermal node N.sub.i in the centroid of mesh element i, a thermal node N.sub.j in the centroid of mesh element j, and an interface thermal node N.sub.ij collocated in the center of the contact interface.

[0018] FIG. 4 is a diagram showing how in the electro-thermal simulation technique disclosed herein, a thermal resistance value r.sub.i is located between node N.sub.i and node N.sub.ij, and a thermal resistance value r.sub.j is located between node N.sub.j and node N.sub.ij.

[0019] FIG. 5 contains graphs showing the correspondence between electro-thermal simulation results using the technique disclosed herein and using a more computationally intensity prior art technique.

DETAILED DESCRIPTION

[0020] The following disclosure enables a person skilled in the art to make and use the subject matter disclosed herein. The general principles described herein may be applied to embodiments and applications other than those detailed above without departing from the spirit and scope of this disclosure. This disclosure is not intended to be limited to the embodiments shown, but is to be accorded the widest scope consistent with the principles and features disclosed or suggested herein.

[0021] First, a more detailed mathematical background explaining existing techniques for electro-thermal modeling will be given. Then, the challenges in applying those techniques will be described. Finally, a novel electro-thermal modeling technique developed by the Inventors will be described in detail.

[0022] To begin the detailed mathematical background, let Ω.Math..sup.3 be a region describing a device formed from thermally conducting materials with n electrical sources. The thermal behavior of the device can be described by the following non-linear heat diffusion equation, for r∈∩ and t>0:

[00003]$\begin{array}{cc}c\left(r,t\right)\frac{\partial \theta }{\partial t}\left(r,t\right)-\nabla .Math.\left(k\left(r,\theta \left(r,t\right)\right)\nabla \theta \left(r,t\right)\right)={.Math.}_{i=1}^n{g}_i\left(r\right)P_i\left(t\right)& \left(1\right)\end{array}$

where θ(r, t) is the unknown temperature rise, c(r, t) is the thermal capacity, k(r, θ(r, t)) is the thermal conductivity, {P.sub.i(t)} are the input powers injected into each source, {g.sub.i(r)} are shape functions describing how powers are distributed in Ω. Observe that the problem is non-linear since the thermal conductivity k(r, θ(r, t)) is temperature-dependent. Following the approach described in document [11], assuming the non-linearity to have the form of k(r, θ(r, t))=k.sub.0(r)λ(r, t)+k.sub.0(r), where the auxiliary variable λ(r,t) is expressed as:

[00004]$\begin{array}{cc}\lambda \left(r,t\right)={\left(1+\frac{\theta \left(r,t\right)}{T_0}\right)}^{\sigma \left(r\right)}-1& \left(2\right)\end{array}$

[0023] where k.sub.0(r) and σ(r) are material dependent parameters and T.sub.0 is the ambient temperature. Adopting this non-linearity in equation (1) yields:

[00005]$\begin{array}{cc}c\left(r,t\right)\frac{\partial \theta }{\partial t}-\nabla .Math.\left[\left(k_0\left(r\right)\lambda +k_0\left(r\right)\right)\nabla \theta \right]={.Math.}_{i=1}^n{g}_i\left(r\right)P_i\left(t\right)& \left(3\right)\end{array}$$\begin{array}{cc}\left(T_0+\theta \right)\frac{d\lambda }{\mathrm{dt}}=\sigma \left(r\right)\left(1+\lambda \right)\frac{d\theta }{\mathrm{dt}}& \left(4\right)\end{array}$

where the unknowns are θ=θ(r, t) and λ=λ(r, t).

[0024] The equation can then be discretized in space by considering a mesh of m elements with a thermal node collocated in each centroid. Several methods may be utilized for this purpose, such as Cell-Center Finite Volumes (see, document [15]). After the discretization, the equation is transformed into a large matrix, as m easily exceeds 10.sup.6. The direct solution of such a system is too processing intensive for numerical solvers (see, document [2]), and therefore an accurate Dynamic Compact Thermal Model (DCTM) is to be derived.

[0025] Here, a Model Order Reduction (MOR) approach (see, document [11]) based on Multi-Point Moment Matching is utilized, where the full system to be solved is projected onto the {circumflex over (m)}-dimensional Krylov subspace spanned by the first terms of a Volterra's series expansion of the Laplace transform of the solutions θ(t) and λ(t) through projection matrices V, W∈.sup.m×{circumflex over (m)}, with {circumflex over (m)}<<m.

[0026] A reduced system is thus obtained through Galerkin projection and changing the variables as such:

θ(t)=Vθ(t) λ(t)=W{circumflex over (λ)}(t)   (5)

[0027] The reduced system is then solved by an iterative alternating procedure, fixing one of the unknowns in each of the two equations and solving them with respect to each other during each iteration. Once the solution to the reduced system has been found, original solutions to θ(t) and λ(t) are recovered through equation (5).

[0028] When applying the above technique, however, two challenges are to be addressed. The first challenge concerns the discretization of the non-linear terms of equations (3) and (4), namely

[00006]$\nabla .Math.\left(k_0\left(r\right)\lambda \nabla \theta \right),\theta \frac{\partial \lambda }{\partial t}\mathrm{and}\sigma \left(r\right)\lambda \frac{\partial \theta }{\partial t}.$

[0029] When Ω is non-homogenous, adjacent mesh elements may correspond to different materials that have different thermal behaviors (for example, silicon and oxide). In this case, material dependent parameters k.sub.0(r) and σ(r) vary from one mesh element to another, and therefore the unknowns θ and λ vary as well. To properly recover the non-linear terms, information about the temperature at the contact surface between different elements is to be known, which is not the case when utilizing standard Finite Volume discretization. To handle this issue, a specific discretization procedure able to accurately describe at the same time the non-linear terms and the temperature rises on interfaces between materials is needed, and such discretization procedure will be described below.

[0030] Before, however, the second challenge is described. This challenge arises when solving the DCTM through the procedure described above. When {circumflex over (λ)} is fixed to solve the system with respect to {circumflex over (θ)}, the system dimensions can be reduced to {circumflex over (m)} using MOR techniques. However, the non-linear term is still expressed in terms of the full λ, as will be explained below. This makes the resolution of even the reduced system to be quite computationally intensive. Therefore, a further reduction of the system that focuses only on the non-linear terms is needed.

[0031] The novel discretization procedure for the non-linear terms in electro-thermo modeling of devices containing non-homogenous materials, developed by the Inventors, is now described.

[0032] First, the hardware 10 used to perform the electro-thermo modeling of this disclosure will be described with reference to FIG. 1. The hardware 10 may be a computing device, and includes a microprocessor 11 that reads instructions and data from, and writes instructions and data to, a volatile memory 12 and a non-volatile memory 13. The computing device 10 includes a display 15 used by the microprocessor 11 to display data to the user. A network interface 14 is used by the microprocessor 11 to retrieve data from, and send data to, a local area network and/or a wide area network. The computing device 10 includes user interface devices 16 used by the microprocessor 11 to obtain user input.

[0033] The non-volatile memory 13 and volatile memory 12 cooperate with the microprocessor 11 to perform the electro-thermo modeling techniques described hereinbelow, with the microprocessor 11 executing instructions stored in the non-volatile memory 13 and/or the volatile memory 12.

[0034] With reference to FIG. 2, in general, first input parameters about the device to be modeled are accepted (Block 21), and then a novel finite volume discretization of the non-linear heat diffusion problem by adding novel interface nodes to the traditional set of inner nodes is performed to reformulate the non-linear heat diffusion problem (Block 22), enabling the electro-thermo modeling of devices containing non-homogenous materials. The non-linear term of the problem is reformulated as the application of a smaller operator, depending upon an auxiliary variable, which serves to reduce memory and computational cost of the technique (Block 23). Then, a MOR technique is combined with Hyper Reduction of the reformulated nonlinear term (Block 24) to enable recovery of the full solution (Block 25). It is to be understood that all steps and calculations described below are performed by cooperation between the various components of the hardware 10.

[0035] A. Finite Volume Discretization

[0036] First, consider the heat diffusion problem in the presence of non-linear thermal conductivity behavior, particularly when the thermal conductivity is not constant and is dependent upon temperature. The general form of the non-linear heat diffusion problem to solve is:

[00007]$c\left(r,t\right)\frac{\partial \theta }{\partial t}\left(r,t\right)-\mathrm{div}\left(k\left(r,\theta \left(r,t\right)\right)\nabla \theta \left(r,t\right)\right)=g\left(r,t\right)-k\left(r,\theta \left(r,t\right)\right)\frac{\partial \theta }{\partial n}\left(r,t\right)=h\left(r,\theta \left(r,t\right)\right)\theta \left(r,t\right)$

where θ(r, t) is the unknown temperature rises in the device being modeled, c(r, t) is the thermal capacity, k(r, θ, r, t)) is the non-linear thermal conductivity, h is the boundary heat exchange coefficient, and g(r) is the shape function modeling the electrical behavior of the device being modeled.

[0037] In document [11], the following auxiliary variable is introduced:

[00008]$\lambda \left(r,t\right)={\left(1+\frac{\theta \left(r,t\right)}{T_0}\right)}^{\sigma \left(r\right)}-1$

[0038] The introduction of λ(r,t) allows rewriting of the heat diffusion problem so that only quadratic non-linearities occur:

[00009]$c\left(r,t\right)\frac{\partial \theta }{\partial t}\left(r,t\right)-\mathrm{div}\left[\left(k_0\left(r\right)\lambda \left(r,t\right)+k_0\left(r\right)\right)\nabla \theta \left(r,t\right)\right]=\underset{i=1}{\overset{n}{.Math.}}{g}_i\left(r\right)P_i\left(t\right)\left(T_0+\theta \left(r,t\right)\right)\frac{\partial \lambda }{\partial t}\left(r,t\right)=\sigma \left(r\right)\left(1+\lambda \left(r,t\right)\right)\frac{\partial \theta }{\partial t}\left(r,t\right)$

[0039] The unknowns then become θ(r,t) and λ(r,t), and the thermal capacity parameter is c(r,t) while the material dependent thermal conductivity parameter is k.sub.0(r).

[0040] The problem can be discretized as follows. The domain of the device to be modeled can be split into m tetrahedrons. An inner thermal node N.sub.i can be collocated in the centroid of each mesh element (tetrahedron) with i=1, . . . ,m. Then, an additional interface node N.sub.ij can be collocated in the center of each contact surface between an N.sub.i tetrahedron and an N.sub.j tetrahedron. This is visually represented in FIG. 3.

[0041] Then, the interconnection structure is derived from the coupled electrical problem, as shown in FIG. 4. For each additional interface node N.sub.ij, two more thermal resistance values are considered, namely the resistance r.sub.i occurring in the connection between the inner node N.sub.i and the interface node N.sub.ij, and the resistance r.sub.j occurring in the connection between the interface node N.sub.ij and the inner node N.sub.j. Overall therefore, this results in a total number of inner thermal nodes plus a total number of interface thermal nodes, this sum being denoted as m.sub.e.

[0042] The discretized non-linear heat diffusion problem then becomes:

[00010]$C\frac{\partial }{\partial t}\theta \left(t\right)+K\theta \left(t\right)+\Delta K\left(\theta \left(t\right).Math.\lambda \left(t\right)\right)=\mathrm{GP}\left(t\right)T_{0}N\frac{\partial \lambda }{\partial t}\left(t\right)+\Delta N\left(\theta \left(t\right).Math.\frac{\partial \lambda }{\partial t}\right)=M\frac{\partial \theta }{\partial t}\left(t\right)+\Delta M\left(\frac{\partial \theta }{\partial t}\left(t\right).Math.\lambda \left(t\right)\right)$

[0043] The matrices C, K, and G respectively describe thermal capacity, thermal conductivity between each mesh element, and the fraction of power dissipated in each node. P(r) is the vector containing the electrical powers for each source.

[0044] The matrix C is computed as a diagonal matrix of thermal capacity, with entry c.sub.ii≠0 if the index i refers to an inner node, with i=1, . . . ,m.sub.e.

[0045] The matrix G is computed as a diagonal matrix of dissipated power, with entry g.sub.ii≠0 if the index i refers to an inner node, with i=1, . . . ,m.sub.e.

[0046] The matrix K is computed as a m.sub.e×m.sub.e matrix, with k.sub.ij=−r.sub.i if i≠j and the resistance occuring in the connection is between N.sub.i and N.sub.ij, and k.sub.ij=the sum in absolute value of the elements of the row if i=j.

[0047] The matrix N is an m×m identity matrix, and M is a diagonal m×m matrix where m.sub.ij=σ(N.sub.i), with σ being a material dependent function.

[0048] B. Reformulation of the Non-Linear Term of the Problem

[0049] The non-linear term in the now discretized problem is given by:

ΔK(θ(t).Math.λ(t))

with λ(r,t) being fixed and independent from r, or, equivalently, λ(r,t)=λ(t).

[0050] The nonlinear matrix ΔK has a size of m.sub.e×(m.Math.m.sub.e) and is sufficiently large such that it is not desirable to be stored in the volatile memory or employed in computations. As such, the nonlinear term is rewritten by finding a squared matrix A.sub.λ(t) of size m.sub.e×m.sub.e such that:

ΔK(θ(t).Math.λ(t))=A.sub.λ(t)θ(t)

A.sub.λ(t) is computed as follows. Each diagonal value k.sub.ii of the matrix K is replaced with the non-linear resistance value λ(N.sub.i,t)k.sub.ii. Each term λ(N.sub.i,t) is computed in the inner nodes as follows:

[00011]$\lambda \left(N_i,t\right)={\left(1+\frac{\theta \left(N_i,t\right)}{T_0}\right)}^{\sigma \left(N_i\right)}-1$

[0051] It can be observed that A.sub.λ(t) depends on the discretization described herein and on the thermal conductivity matrix K. This reformulation of A.sub.λ(t) allows simplification of the non-linear term as the operator A.sub.λ(t) has a size of m.sub.e×m.sub.e while ΔK has a size of m.sub.e×(m.Math.m.sub.e)—this step is particularly useful and novel, as it greatly reduces the computing and memory overhead used by the simulation techniques described herein, and greatly reduces the computation time for such simulation techniques.

[0052] C. Hyper Reduction of the Reformulated Non-Linear Term

[0053] Non-linear MOR techniques are dimensionality reduction techniques that compute projection matrices V,W able to derive a compact form for the discretized non-linear heat diffusion problem:

[00012]$\hat{C}\frac{\partial }{\partial t}\hat{\theta }\left(t\right)+\hat{K}\hat{\theta }\left(t\right)+\Delta \widehat{K}\left(\hat{\theta }\left(t\right).Math.\hat{\lambda }\left(t\right)\right)=\hat{G}P\left(t\right)T_0\hat{N}\frac{\partial \hat{\lambda }}{\partial t}\left(t\right)+\Delta \hat{N}\left(\hat{\theta }\left(t\right).Math.\frac{\partial \hat{\lambda }}{\partial t}\right)=\hat{M}\frac{\partial \hat{\theta }}{\partial t}\left(t\right)+\Delta \hat{M}\left(\frac{\partial \hat{\theta }}{\partial t}.Math.\hat{\lambda }\left(t\right)\right)$

[0054] The matrix V has a size of {circumflex over (m)}×m.sub.e and W has a size of {circumflex over (m)}×m, where m.sub.e>m>>{circumflex over (m)}. Reduced matrices are then defined as follows:

Ĉ=V.sup.TCV {circumflex over (m)}×{circumflex over (m)}

{circumflex over (K)}=V.sup.TKT {circumflex over (m)}×{circumflex over (m)}

Δ{circumflex over (K)}=V.sup.TΔK(V.Math.W) {circumflex over (m)}×{circumflex over (m)}.sup.2

Ĝ=V.sup.TG {circumflex over (m)}×{circumflex over (m)}

{circumflex over (N)}=W.sup.TNW {circumflex over (m)}×{circumflex over (m)}

{circumflex over (M)}=W.sup.TMV {circumflex over (m)}×{circumflex over (m)}

Δ{circumflex over (N)}=W.sup.TΔN(V.Math.W) {circumflex over (m)}×{circumflex over (m)}.sup.2

Δ{circumflex over (M)}=W.sup.TΔM(V.Math.W) {circumflex over (m)}×{circumflex over (m)}.sup.2

[0055] The calculation of the matrices ΔK, ΔM, and ΔN is as follows:

ΔK=Σ.sub.i=1.sup.mA.sub.e.sub.i.Math.e.sub.i.sup.T

ΔM=Σ.sub.i=1.sup.mM.Math.e.sub.i.sup.T

ΔN=Σ.sub.i=1.sup.mN.Math.e.sub.i.sup.T

where e.sub.i=(0, . . . , 0, 1, 0, . . . , 0).sup.T is the ith vector of the basis of the m-dimensional Euclidean space.

[0056] The projection matrices V and W are computed using a non-linear MOR technique based upon Moment Matching and Krylov subspaces. Here, the matrices V and W span a set of solutions of the non-linear problem evaluated in a set of frequencies in the Laplace domain.

[0057] The reduced non-linear matrix ΔK has a size of {circumflex over (m)}×{circumflex over (m)}.sup.2, and Hyper Reduction techniques are employed to further reduce the number of columns in the reduced non-linear matrix ΔK. Through this, a subset of the original full mesh volumes is identified and appropriate weights are applied such that the weighted reduced mesh volumes well approximate the original problem.

[0058] In particular, Energy-Conserving Sampling and Weighting (ECSW) is applied to the reduced non-linear matrix ΔK. It begins with a collection of s snapshots of the solution to the non-linear problem in a set of complex frequencies in the Laplace domain, resulting in a matrix U having dimensions of m.sub.e×s and a matrix A having dimensions of m×s.

[0059] Through the definition of ΔK, construction of a matrix H representing the unassembled contributions of the non-linear forces and a vector b occurs as follows:

H.sub.ie=V.sup.T*∧(e,i)*A.sub.pe*U(:,i),

i=1, . . . , s

e=1, . . . , m

b=H*(1, . . . , 1).sup.T

[0060] Resolution of the following sparse minimization problem with a non-negative variant of the orthogonal matching pursuit algorithm is then performed as:

ξ*=argmin.sub.ξ∥Hξ−b∥.sub.F s.t.ξ≥0 and ∥ξ∥.sub.0<w

[0061] The solution of the problem ξ* represents the sparse vector of weights, and a set E can then be defined as the indices of the non-zero entries of ξ*. Computation of the non-linear model order reduced term is performed by evaluating it on the volumes indexed by E and weighted with relative ξ*, as follows:

Δ{circumflex over (K)}({circumflex over (θ)}(t).Math.{circumflex over (λ)}(t))≈(Σ.sub.e=1.sup.wξ*(e)*W(e,:)*{circumflex over (λ)}(t)*V.sup.TA.sub.pe*V)*{circumflex over (θ)}(t)

[0062] To further increase the computation speed, a singular value decomposition (SVD) is performed on the terms V.sup.T*A.sub.pe*V for e=1, . . . , w, further reducing the complexity of the computations.

[0063] To recover the full solution (i.e., the equivalent simulation result as if the full prior art simulation technique was used), the Hyper Reduced problem is solved and the reduced solution is computed: {circumflex over (θ)}(t) and {circumflex over (λ)}(t). Then, the full solution is recovered from the reduced solution, using the projection matrices V and W:

θ(t)=V{circumflex over (θ)}(t) λ(t)=W{circumflex over (λ)}(t)

[0064] A comparison between the results using the DCTM technique described herein and a prior art finite volume model was made. Using the DCTM technique, 65 minutes of computing time was used to compute the projection matrices V, W, and the vector ξ*, and this calculation is performed but once since it depends solely on the domain Ω. Multiple thermal simulations for different injected powers of Pi were performed, each taking but 12.160 seconds to perform. For comparison, the prior art finite volume model required 197 minutes of computing time. An accuracy comparison between the DCTM technique described herein and the prior art finite volume model can be seen in FIG. 5, where it can be observed that the results are nearly identical.

[0065] D. Improvements to Electro-Thermo Simulation Technology and Computing Devices Performing Electro-Thermo Simulations

[0066] In greater detail, the structure of Smart Power devices is, to date, primarily based on silicon. However, silicon is not the only material utilized in such devices—copper and aluminum are employed in the interconnections of the device, while different oxides are used to isolate different regions. In particular, silicon oxide is adopted to create microstructures within the silicon parts, such as Deep Trench Isolations (DTI) and Silicon on Insulator (SOI), to improve the electrical properties of the device.

[0067] On the other hand, silicon oxide is a poor heat conductor. As a consequence, the adoption of such solutions and the significant increase of the power density (due to the miniaturization of the devices) result in a drastic rise of the operation temperatures and in the formation of high temperature gradients within the devices.

[0068] Therefore, during the design stages of these devices, thermal and electrothermal analysis resulting in an accurate description of the temperature rises within the device are strongly desired to help guarantee its correct operation.

[0069] As explained earlier, the structure of these devices, formed of different materials each with a different thermal behavior, results in a thermal problem which is non-linear and difficult to model with standard prior art 3D-discretization methods (e.g., Cell-Centered Finite Volumes). In fact, such techniques fail when describing non-homogeneous structures with significant non-linearities in their materials, thereby not allowing an accurate analysis of sudden temperature variations within the discretized domain.

[0070] Conversely, the DCTM techniques described herein provide an extremely accurate description of the non-linear behavior of the materials forming the device, regardless of their structure. This allows the realization of accurate simulations even under the high temperatures at which Smart Power devices work, which may not be possible with prior art techniques. Moreover, due to the application of Model Order Reduction (MOR), the extraction of dynamic compact thermal models has been accomplished, resulting in simulations which are not only accurate but also extremely fast, as also explained earlier.

[0071] The Smart Power industry is currently facing significant changes and challenges—new materials like gallium nitride (GaN) and silicon carbide (SiC) will be integrated in Smart Power devices, making the resulting structure extremely complex to analyze and the temperature rise reached during operation even higher.

[0072] Therefore, the discretization techniques described herein not only drastically improve over prior art simulation frameworks, but provide for adequate accuracy to handle the analysis on future devices.

[0073] To better understand the usefulness of the accurate electro-thermal analysis provided by the techniques described herein, consider the simple example of a short-circuit control system of a power output of a voltage regulator. The system needs to be able to detect a short-circuit downstream of the output. Therefore, the output current has to be assessed and the power supply stopped in order to avoid damages to the device. When a short-circuit occurs, the current used in the device significantly increases and, as a consequence, so does the dissipated power.

[0074] The sizing of the device is aimed at guaranteeing its functioning during its operational life. A sudden increase of the dissipated power can force the device to a very high temperature rise for a short time span (before the interruption of the power supply). Such temperature rise may come close to the critical temperature of the device. The critical temperature is a threshold which, when is overcome (due to the leakage effects of the current), results in a further increase of the dissipated power. This, due to the exponential form of the leakage currents, can cause destructive phenomena. The value of the critical temperature is different for each device, but is known a priori.

[0075] Electro-thermal simulations allow forecasting of the behavior of the device, however, it is important to provide an accurate description of the thermal model in order to obtain reliable results.

[0076] Given the exponential form of the electrical effects at stake, if a model overestimates the temperature rises during the simulation, the resulting device will be oversized, employing an unnecessary area, and increasing the final cost of the product. On the other hand, if a model underestimates the temperature rises, the device will be undersized, and therefore more subject to be damaged when a short-circuit occurs during its functioning.

[0077] Both scenarios are to be prevented: in one case the device cost would be too high, hindering competitiveness on the market, while in the other the expenses to redesign the device and the management of the return from the field costs would cause a significant economical loss and image damage for the company. The electro-thermal simulation and analysis techniques described herein facilitate the prevention of both scenarios, using fewer computing and memory resources, in a shorter period of time, representing an advance in electro-thermal modeling technology itself, as well as representing an improvement in the operation of workstation computers performing electro-thermal modeling.

[0078] The introduction of additional thermal nodes on the surfaces of contact between mesh elements in the DCTM techniques described herein leads to several advantages. First of all, it allows description of the non-linear behavior of different materials (e.g., silicon and oxide). Although introducing additional nodes increases the dimension of the problem, the resulting structure of the matrices defining heat diffusion is computationally easier to handle. In fact, the new matrices present a large amount of non-zero elements and their structure allows linearization of the heat diffusion equation, drastically reducing the complexity of the algorithm in both the MOR and the Hyper Reduction stages.

[0079] Listing of documents referenced herein:

[0080] [1] Vladimir Szekely, “On the representation of infinite-length distributed rc one-ports,” IEEE Transactions on Circuits and Systems, vol. 38, no.7, 1991.

[0081] [2] Olek C Zienkiewicz, Robert L Taylor, and Jian Z Zhu, “The finite element method: its basis and fundamentals,” Elsevier, 2005.

[0082] [3] Wilhelmus H A Schilders, Henk A Van der Vorst, and Joost Rommes, “Model order reduction: theory, research aspects and applications,” Springer, 2008.

[0083] [4] Anindya Chatterjee, “An introduction to the proper orthogonal decomposition,” Current Science, vol. 78, no. 7, 2000.

[0084] [5] Francisco Chinesta, Pierre Ladeveze, and Elias Cueto, “A short review on model order reduction based on proper generalized decomposition,” Archives of Computational Methods in Engineering, vol. 18, no. 4, 2011.

[0085] [6] Pieter Jacob Heres, “Robust and efficient Krylov subspace methods for model order reduction,” Ph.D. thesis, Technische Universiteit Eindhoven, 2005.

[0086] [7] Yehea I Ismail, “Efficient model order reduction via multi-point moment matching,” 2004, U.S. Pat. No. 6,789,237.

[0087] [8] Lorenzo Codecasa, Dario D'Amore, and Paolo Maffezzoni, “Compact modeling of electrical devices for electrothermal analysis,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 50, no. 4, 2003.

[0088] [9] Lorenzo Codecasa, Vincenzo d′Alessandro, Alessandro Magnani, Niccolo Rinaldi, and Peter J Zampardi, “Fast novel thermal analysis simulation tool for integrated circuits (fantastic),” in International Workshop on Thermal Investigations of ICs and Systems, 2014.

[0089] [10] Lorenzo Codecasa, Dario D′Amore, and Paolo Maffezzoni, “An arnoldi based thermal network reduction method for electro-thermal analysis,” IEEE Transactions on Components and Packaging Technologies, vol. 26, no. 1, 2003.

[0090] [11] Lorenzo Codecasa, Vincenzo d′Alessandro, Alessandro Magnani, and Niccolo Rinaldi, “Fast nonlinear dynamic compact thermal modeling with multiple heat sources in ultra-thin chip stacking technology,” IEEE Transactions on Components, Packaging and Manufacturing Technology, vol. 7, no. 1, 2016.

[0091] [12] Lorenzo Codecasa, Vincenzo d′Alessandro, Alessandro Magnani, and Niccolo Rinaldi, “Novel approach for the extraction of nonlinear compact thermal models,” in International Workshop on Thermal Investigations of ICs and Systems, 2017.

[0092] [13] RE Mortensen, “Nonlinear system theory: The volterra/wiener approach,” SIAM Review, vol. 25, no. 3, 1983.

[0093] [14] Charbel Farhat, Philip Avery, Todd Chapman, and Julien Cortial, “Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency,” International Journal for Numerical Methods in Engineering, vol. 98, no. 9, 2014.

[0094] [15] Robert Eymard, Thierry Gallouet, and Raphaele Herbin, “Finite volume methods,” in Handbook of numerical analysis, vol. 7. Elsevier, 2000.

[0095] [16] Edoardo Amaldi, Viggo Kann, et al., “On the approximability of minimizing nonzero variables or unsatisfied relations in linear systems,” Theoretical Computer Science, vol. 209, no. 1, 1998.

[0096] [17] Yagyensh Chandra Pati, Ramin Rezaiifar, and Perinkulam Sambamurthy Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” Asilomar conference on signals, systems and computers, 1993.

[0097] It is to be noted that documents [1] through [17] are each incorporated by reference in their entirety.

[0098] Further details may be found in U.S. Pat. No. 9,384,315, entitled “Method, system and computer program product for electrical and thermal analysis at a substrate level”, issued on Jul. 5, 2016, the contents of which are incorporated by reference in their entirety.

[0099] While the disclosure has been described with respect to a limited number of embodiments, those skilled in the art, having benefit of this disclosure, will appreciate that other embodiments can be envisioned that do not depart from the scope of the disclosure as disclosed herein. Accordingly, the scope of the disclosure shall be limited only by the attached claims.