Chip thermal management method based on precomputation random steady-state thermal analysis solver

Abstract

The invention belongs to the technical field of integrated circuit computer aided design / electronic design automation, and relates to a chip active thermal management method adopting a rapid and accurate temperature predictor. According to the method, a steady-state thermal analysis solver based on random walk suitable for a mixed boundary condition is improved by utilizing the characteristics that only a power diagram at a heat source changes and a geometrical shape, material parameters and boundary conditions are unchanged in a thermal management process, and a changed power item is separated from other fixed parts; in the off-line stage, fixed geometrical shapes, material parameters and boundary conditions are pre-calculated and stored, and microsecond-level rapid prediction of the hot-spot temperature is achieved through vector multiplication during chip operation. According to the method, partial inversion is carried out on a large-scale matrix constructed by a compact thermal model based on finite difference by utilizing a random method. And task mapping and dynamic voltage frequency regulation measures are combined, so that extremely high chip working frequency can be kept while all thermal conflicts are eliminated. The method has the advantages of no discretization, high calculation efficiency, natural parallelism and the like.

Description

Technical Field

[0001] This invention belongs to the technical field of integrated circuit computer-aided design / electronic design automatic CAD / EDA, specifically relating to a proactive thermal management method for chips based on a pre-computation stochastic steady-state thermal analysis solver. This method can maintain a high chip operating frequency while eliminating all thermal conflicts; simultaneously, it has advantages such as no need for discretization, high computational efficiency, and high parallel efficiency in the offline pre-computation stage. Background Technology

[0002] As the feature size of semiconductor manufacturing processes shrinks to the nanometer level, technologies such as three-dimensional integrated circuits (3D-IC) and chiplets have developed rapidly to meet the demand for continuously improving performance. However, the rapidly increasing number of transistors causes the power density and temperature of the chip to rise rapidly, leading to serious thermal problems that affect the chip's performance, reliability, and aging [1]. Therefore, heat dissipation has become one of the key challenges in the design and manufacturing of chips such as 3D-IC and chiplets. Achieving accurate and fast integrated circuit thermal simulation is of great significance in the early stages of chip design and chip operation.

[0003] Thermal analysis or thermal simulation is an important sub-problem in integrated circuit thermal problems, focusing on calculating the temperature field distribution by solving the heat conduction equation. Currently, deterministic methods involving full-space discretization and solving large-scale linear equations are widely used in chip thermal analysis problems. Deterministic methods mainly include finite-element methods (FEMs) and finite-difference methods (FDMs). Finite-element methods, such as the widely used ANSYS[2] and COMSOL[3] thermal simulators in industry, discretize the entire region and then construct and solve the linear equation system. They can achieve extremely high accuracy through extremely fine discretization and are very flexible in dealing with complex structures. Finite-difference methods, such as the popular HotSpot[4] and 3D-ICE[5] in academia, use compact thermal models (CTMs) to discretize the region into an equivalent resistor-capacitor (RC) network and solve the RC equation. In handling complex geometries and boundary conditions, they are less flexible than the finite element method. While these deterministic methods can achieve high-precision results, their runtime and memory costs are high, making them unsuitable for the early stages of IC design where repeated thermal simulations and optimizations are required, as well as for real-time temperature prediction during chip operation. Furthermore, they can only obtain local hotspot temperatures by solving the global temperature field, and cannot obtain only local hotspot temperatures, thus being inefficient in acquiring local hotspot temperatures.

[0004] Unlike deterministic methods, stochastic methods based on the Feynman-Kac formula [6] can obtain local hot spot temperatures efficiently and accurately through high parallelism. In stochastic methods, the accurate handling of three boundary conditions, namely Dirichlet (first type), Neumann (second type), and Robin (third type), has an important impact on the final accurate solution of the thermal equation. Dirichlet boundary has a fixed surface temperature and can be easily handled by the path absorption method [7,8]. Neumann boundary is usually empirically regarded as adiabatic boundary in the engineering practice of thermal problems and is also easy to handle. Robin boundary, also known as convection boundary, is one of the key ways for the chip to exchange heat with the outside, and is difficult to handle in stochastic methods. In recent years, reference [9] proposed a random walk method combining Walk-On-Grids (WOG) and Walk-On-Cuboids (WOC) to solve steady-state thermal analysis problems. The Neumann boundary is handled by the reflection method commonly used in engineering, and the Robin boundary is handled by an empirical heuristic algorithm that sets virtual Dirichlet nodes and probabilistic absorption or reflection. However, the calculated temperature field has a large error, with an error of more than 20℃ in some cases. Reference

[10] proposed a Path Integral Random Walk (PIRW) solver, which combines an accurate local time allocation scheme

[11] and the Feynman-Kac functional.

[12] It accurately handles the Neumann and Robin boundary conditions that are prevalent in thermal problems. However, although the PIRW solver achieves a speedup of more than 100 times compared to the finite element method, its simulation speed of seconds is still far from meeting the needs of real-time temperature prediction of chips.

[0005] Thermal management focuses on real-time temperature control during chip operation. Passive methods detect thermal conflicts through temperature sensors on the chip and take measures such as dynamic voltage frequency adjustment, task mapping, and task migration

[13] to reduce the chip operating temperature. In contrast, active methods predict the future temperature of the chip and take the above preventive measures before thermal conflicts occur. When the chip is running under high load, active thermal management will be triggered periodically to predict the future temperature of the chip hotspots to determine whether cooling measures are needed. The effectiveness of active thermal management depends heavily on the accuracy of its internal temperature predictor, which also needs to be very efficient and lightweight so that it can be frequently called in real time.

[0006] Because active thermal management operates in real time, it requires extremely high execution speeds from the built-in temperature predictor, typically on the order of milliseconds or even microseconds. Therefore, deterministic thermal analysis methods such as the finite element method cannot be simply used. Instead, data-driven empirical models or modified physical model methods are usually employed.

[0007] Data-driven empirical models include autoregressive moving average models

[14] , neural networks

[13] , and regression models

[15] . In 2009, Coskun et al. proposed an autoregressive moving average (ARMA) model to predict future temperatures. The ARMA model is modeled as a stationary stochastic process, including an offline setup process and an online real-time adjustment

[14] . In 2010, Ge et al. proposed a neural network method to predict future temperatures, with an average performance 79% better than the ARMA model and a maximum predicted temperature error of 2.5℃

[13] . In 2019, et al. proposed a temperature prediction method for a regression model. In the offline stage, the power and temperature of the actual operation are used for training. In the online prediction stage, it introduces an error correction mechanism for looking up and updating the error table to improve the online prediction accuracy, which reduces the average error of temperature prediction to 1.15℃

[15] .

[0008] The modified physical model is commonly used when the chip structure and physical parameters are known, and a temperature predictor based on the compact thermal model (CTM) is used [4,16]. The compact thermal model is essentially a finite difference method, which uses the chip structure and thermal parameters to construct an equivalent RC circuit and obtains the RC circuit results using standard circuit simulation. When the chip is running, the general geometry and boundary conditions are fixed, and only the power consumption diagram of the heat source area on the right side changes. In essence, the temperature predictor of the compact thermal model CTM is essentially a large-scale matrix inversion offline and a simple matrix-vector multiplication online. When used as a temperature predictor in thermal management methods, CTM is often further simplified, such as placing only one hot node that needs to be measured at the center of each CPU core [17,18,19], approximating the temperature of so-called unmeasurable nodes in non-heat source areas

[17] , introducing empirical parameters to further simplify the thermal network

[18] , ignoring convection heat transfer at the boundary

[19] , etc.

[0009] Based on the current state and foundation of existing technologies, the inventors of this application intend to modify the thermal analysis algorithm based on the stochastic method, and provide a high-precision stochastic method-based temperature predictor pPIRW for active thermal management of chips. This method can serve as a powerful tool for online rapid temperature prediction and management, and has a significant speed advantage after applying large-scale parallel computing technology. This method can be extended to fast field analysis and field prediction problems of Laplace and Poisson equations with mixed boundary conditions in other engineering fields.

[0010] References related to this invention:

[0011] [1]Sheng-Chih Lin and Kaustav Banerjee.2008.Thermal challenges of 3DICs.In Wafer Level 3-D ICs Process Technology.Springer,1–26.

[0012] [2] Inc.ANSYS, "Ansys," https: / / www.ansys.com / , 2018.

[0013] [3] Inc.COMSOL, "Comsol," https: / / www.comsol.com / , 2022.

[0014] [4]K.Skadron, MRStan, K.Sankaranarayanan, W.Huang, S.Velusamy, andD.Tarjan, "Temperature-aware microarchitecture: Modeling and implementation," ACM Trans.Archit.Code Optim., vol.1, no.1, pp.94–125, 2004.

[0015] [5] F.Terraneo, A.Leva, W.Fornaciari, M.Zapater, and D.Atienza, "3D-ICE3.0: efficient nonlinear MPSoC thermal simulation with pluggable heat sinkmodels," IEEE Trans.on CAD, vol.41, no.4, pp.1062–1075, 2021.

[0016] [6]J.-P.Morillon,“Numerical solutions of linear mixed boundary valueproblems using stochastic representations,”International journal fornumerical methods in engineering,vol.40,no.3,pp.387–405,1997.

[0017] [7]Abdolhossein Haji-Sheikh and Ephraim M Sparrow.The floating randomwalk and its application to monte carlo solutions of heat equations.SIAMJournal on AppliedMathematics,14(2):370–389,1966.

[0018] [8]S Talebi,K Gharehbash,and HR Jalali.Study on random walk and itsapplication to solution of heat conduction equation by monte carlomethod.Progress in Nuclear Energy,96:18–35,2017.

[0019] [9]Y.Liang,W.Yu,and H.Qian,“A hybrid random walk algorithm for 3-Dthermal analysis ofintegrated circuits,”inASP-DAC.IEEE,2014,pp.849–854.

[0020]

[10] L.Yang,C.Ding,C.Yan,D.Zhou,and X.Zeng,“A high-precisionstochastic solver for steady-state thermal analysis with Fourier heattransfer robin boundary conditions,”in 202241stACM / IEEE ICCAD,2022.

[0021]

[11] C.Ding,Y.Zhou,W.Cai,X.Zeng,and C.Yan,“A path integral Monte Carlo(PIMC)method based on Feynman-Kac formula for electrical impedancetomography,”Journal ofComputational Physics,p.111862,2023.

[0022]

[12] V.G.Papanicolaou,“The probabilistic solution of the thirdboundary value problem for second order elliptic equations,”ProbabilityTheory and related fields,vol.87,no.1,pp.27–77,1990.

[0023]

[13] Y.Ge,Q.Qiu,and Q.Wu,“A multi-agent framework for thermal awaretask migration in many-core systems,”IEEE Transactions on Very Large ScaleIntegration(VLSI)Systems,vol.20,no.10,pp.1758–1771,2011.

[0024]

[14] A.K.Coskun,T.S.Rosing,and K.C.Gross,“Utilizing predictors forefficient thermal management in multiprocessor SoCs,”IEEE Transactions onComputer-Aided Design of Integrated Circuits and Systems,vol.28,no.10,pp.1503–1516,2009.

[0025]

[15] E.W.W¨achter,C.De Bellefroid,K.R.Basireddy,A.K.Singh,B.M.Al-Hashimi,and G.Merrett,“Predictive thermal management for energy-efficientexecution of concurrent applications on heterogeneous multicores,”IEEETransactions on Very Large Scale Integration(VLSI)Systems,vol.27,no.6,pp.1404–1415,2019.

[0026]

[16] S.Pagani,J.-J.Chen,M.Shafique,and J.Henkel,“MatEx:efficienttransient and peak temperature computation for compact thermal models,”in2015 Design,Automation&Test in Europe Conference&Exhibition(DATE).IEEE,2015,pp.1515–1520.

[0027]

[17] S.Sharifi,R.Ayoub,and T.S.Rosing,“Tempomp:Integrated predictionand management oftemperature in heterogeneous MPSoCs,”in 2012 Design,Automation&Test in Europe Conference&Exhibition(DATE).IEEE,2012,pp.593–598.

[0028]

[18] M.Li,W.Liu,L.Yang,P.Chen,and C.Chen,“Chip temperatureoptimization for dark silicon many-core systems,”IEEE Transactions onComputer-Aided Design ofIntegrated Circuits and Systems,vol.37,no.5,pp.941–953,2017.

[0029]

[19] G.Bhat,G.Singla,A.K.Unver,and U.Y.Ogras,“Algorithmic optimizationof thermal and power management for heterogeneous mobile platforms,”IEEETransactions on Very Large Scale Integration(VLSI)Systems,vol.26,no.3,pp.544–557,2017.

[0030]

[20] Y.Zhou,W.Cai,and E.Hsu,“Computation of the local time ofreflecting brownian motion and the probabilistic representation of theneumann problem,”Communications in Mathematical Sciences,vol.15,no.1,pp.237–259,2017.

[0031]

[21] P.Hsu, Reflecting Brownian motion, boundary local time and the Neumann problem. Stanford University, 1984.

[0032]

[22] Eigen, "Eigen c++template library," http: / / eigen.tuxfamily.org / , 2021. Summary of the Invention

[0033] The purpose of this invention is to propose a chip thermal management method based on a pre-computation stochastic steady-state thermal analysis solver, building upon the current state and foundation of existing technologies. This method aims to achieve a high-precision, high-speed temperature predictor in proactive chip thermal management. The method decomposes the chip thermal management problem into an offline thermal analysis (prediction) problem and an online dynamic thermal management optimization problem.

[0034] This method consists of two stages: offline computation and online thermal management. In the offline computation stage, a steady-state thermal analysis solver based on pre-calculation path integral random walk (pPIRW), suitable for mixed boundary conditions (mixed BCs), is used for pre-computation and storage. During online thermal management at chip runtime, vector multiplication is used to achieve microsecond-level rapid prediction of chip hotspot temperatures, and active thermal management techniques such as task mapping and dynamic voltage and frequency scaling (DVFS) are employed to reduce the chip's maximum temperature. Experimental results show that this method can maintain a high chip operating frequency while eliminating all thermal conflicts; furthermore, the offline pre-computation stage offers advantages such as no discretization required, high computational efficiency, and high parallel efficiency.

[0035] This invention addresses the problem of active thermal management in chips by modifying a stochastic thermal analysis algorithm and proposing a high-precision stochastic temperature predictor, pPIRW, for active thermal management. Since the geometry and boundary conditions of integrated circuits are fixed during chip operation, the only variable is the real-time power graph. Therefore, based on the high-precision PIRW solver, the changing real-time power factor is extracted from the temperature expression, and the remaining relevant terms, namely the decaying Feynman-Kac functionals of the Robin boundary in all paths, are pre-calculated through numerous random walks. This time-consuming process is considered to be performed offline, while the remaining online temperature prediction is achieved through simple vector multiplication, ultimately resulting in a fast real-time temperature prediction at the microsecond level.

[0036] Specifically,

[0037] In this invention:

[0038] In the offline thermal analysis phase, the chip thermal analysis problem is transformed into a Poisson equation solution problem with mixed boundary conditions. Based on the Feynman-Kac formula, the Walk-On-Spheres (WOS) method of standard reflection Brownian motion is used to absorb or reflect the first, second, and third types of boundaries, respectively, to obtain the probabilistic solution of the hot spot temperature in the chip steady-state thermal analysis problem. Furthermore, the power term is decoupled from the random walk path, and pre-calculation is performed through random walks to extract the varying power term from other invariant parameters in the chip thermal management problem.

[0039] Extracting from the portion, we obtain k' hotspots (H1, H2, ..., H...). k' Pre-computed matrix A and vector at position ) And a temperature prediction formula in the form of fast matrix multiplication and addition is obtained. In the online dynamic thermal management stage, the pre-calculated...

[0040] Matrix A and vector It is applied to active thermal management methods to achieve real-time temperature prediction, and uses task mapping and dynamic voltage frequency regulation (DVFS) measures to control the chip hot spot temperature to not exceed the set temperature threshold.

[0041] This invention discloses a chip thermal management method based on a pre-computed path integral (pPIRW) random walk solver, which solves the active thermal management problem of chips. Its flowchart is shown below. Figure 1 As shown.

[0042] Input parameters: region Ω and its boundaries Geometric information, material thermal conductivity parameters, ambient temperature T amb ,boundary The three types of boundaries Γ i Boundary conditions φ i (i = 1, 2, 3), number of Monte Carlo sampling paths N, maximum number of path jumps NP, path termination error tolerance ε, k' hotspot locations H1, H2, ..., H k' D near the boundaries of the second and third categories ò The thickness of the region ò and the jump radius Δx of the random walk inside, the width ò1 of the absorption band near the first-class boundary, and the power diagram of the heat source region over time.

[0043] Output parameters: Predict the steady-state temperature of the chip hotspot and adjust the task mapping and dynamic voltage frequency regulation DVFS level based on the predicted temperature.

[0044] The specific steps include:

[0045] Step 1: Decompose the chip's thermal management problem into a thermal analysis (prediction) problem and an online dynamic thermal management optimization problem. Further, transform the chip thermal analysis problem into a Poisson equation solving problem with mixed boundary conditions.

[0046] Step 2: Based on the Feynman-Kac formula, obtain the probabilistic solution for the hotspot temperature of the chip steady-state thermal analysis problem. Decouple the power term from the random walk path to obtain k' hotspots (H1, H2, ..., H...). k' Temperature prediction formula in matrix multiplication and addition form at ().

[0047] Step 3: Initialize a connection from hotspot H i The starting random walk path, i.e., the vector 'a' initialized in the pre-computation process. and boundary contribution value Record array.

[0048] Step 4: Based on the current point X t At the given location, perform the corresponding random walk procedure, using WOS random walk in the passive region; update near the first type Dirichlet boundary. And terminate the random walk directly; calculate and update the local time near the Neumann boundary of the second kind. Calculate and update the local time near the third type of Robin boundary. and In the active region, a WOG random walk is performed, and the corresponding element of vector a is updated according to the grid position number of the current point.

[0049] Step 5: Determine X after the jump t+1 The position of X t+1 Located on the absorption boundary or If the maximum number of jumps in the path is reached (NP), the random walk ends, and the process proceeds to step 6. Otherwise, the jump to the next path is... t+1 As the current point X t Proceed to step 4.

[0050] Step 6: Determine whether all k' hotspots have completed N random walks. If they have, proceed to step 7; otherwise, proceed to step 3 to continue random walks.

[0051] Step 7: Calculate the sum of vector 'a' for each hotspot after executing N paths. The average value of all hotspots is calculated, and the pre-computed values ​​of all hotspots are assembled into matrix A and vector A.

[0052] Step 8: Combine the pre-computed matrix A and vector It is applied to active thermal management methods for real-time temperature prediction, and uses task mapping and DVFS measures to control the chip hot spot temperature to not exceed the threshold.

[0053] In step 1 of this invention, the thermal management problem of the chip is decomposed into a thermal analysis (prediction) problem and an optimization problem of online dynamic thermal management. Furthermore, the thermal analysis problem of the chip is transformed into a Poisson equation solving problem with mixed boundary conditions.

[0054] The chip thermal analysis problem that this invention aims to solve is as follows: Figure 2 As shown. Figure 2 This demonstrates how chip thermal analysis can be applied to real-world chip packaging. Figure 2 As shown in a, the simplified thermal model (such as...) Figure 2 b), finally to the mathematical model ( Figure 2 The definition process of c). Figure 2 In model a, the chip portion (red) within the chip package is the heat source area where temperature measurement is required. Other layers, such as heat sinks, heat diffusers, thermal interface materials (TIMs), substrates, and printed circuit boards (PCBs), are modeled as non-heat source areas.

[0055] Typically, the primary heat dissipation path of a chip is through the heatsink (top surface), and therefore it is modeled as a convective (Robin) boundary. Secondary heat transfer paths are through the PCB (bottom surface), which can be modeled as isothermal (Dirichlet), adiabatic (Neumann), or convective (Robin) boundaries, depending on the specific application. The four side surfaces, due to minimal heat conduction, are usually modeled as adiabatic (Neumann) boundaries.

[0056] Figure 2 c presents the mathematical model for the three-dimensional steady-state heat conduction problem. The three-dimensional region for steady-state thermal analysis is Ω, and the temperature field T(x) within the region Ω satisfies the Poisson equation. The boundary of the region... The boundary is a mixed boundary, meaning it contains Dirichlet boundary Γ1, Neumann boundary Γ2, and Robin boundary Γ3, or only one or two of them. Given the three types of boundary conditions, the field differential equations and boundary conditions are as follows:

[0057]

[0058] Where T(x) is the temperature within region Ω, p d φ(x) is the heat generation power density, and k(x) is the thermal conductivity of the material. On the first kind Dirichlet boundary Γ1, the temperature T(x) is a known φ1(x); on the second kind Neumann boundary Γ2, the external normal derivative of temperature T at the boundary is... Given φ2(x), in steady-state thermal analysis, it is usually an adiabatic boundary condition, i.e. The Robin boundary, also known as the convection boundary, corresponds to the interface that exchanges heat with the surrounding fluid, such as the surface of a cooling fan and radiator. The parameter c(x) in the Robin boundary condition is:

[0059]

[0060] Where h(x) is the heat convection coefficient at the Robin boundary, which is typically used in steady-state thermal problems. Where T amb The ambient temperature is known.

[0061] The active thermal management process of a chip can be described as an optimization problem:

[0062] Maximize average V / flevel, (3)

[0063] stT≤T threshold

[0064] The optimization objective is to maximize speed when executing a set of applications, i.e., to maximize the average voltage or frequency (V / f) of the computing core, while keeping the hotspot temperature T at a threshold T. threshold the following.

[0065] In step 2 of this invention, the probabilistic solution for the hotspot temperature of the chip steady-state thermal analysis problem is obtained based on the Feynman-Kac formula. The power term is decoupled from the random walk path to obtain k' hotspots (H1, H2, ..., H...). k' Temperature prediction formula in matrix multiplication and addition form at ().

[0066] Considering the Poisson equation (1) with mixed boundary conditions, according to the Feynman-Kac formula, the probabilistic solution for the temperature value of point x in region Ω is [6]:

[0067]

[0068] Among them, E x Let X be the mathematical expectation. t This refers to Reflecting Brownian Motion (RBM). It is the moment when the reflected Brownian motion reaches the Dirichlet boundary. It is the position corresponding to that moment; Let D be the Feynman-Kac functional, and L(t) be the local time of the reflection Brownian motion near the Neumann and Robin boundaries [20,21], where D ò It is a strip-shaped area with a width of ò near the boundary. It is in D ò The index function within the region. The integral terms in (4)-(6) represent the contributions of the right-hand side of the Poisson equation to the temperature at the point to be measured under the Dirichlet, Neumann, and Robin boundary conditions, respectively.

[0069] Reflective Brownian motion is usually achieved by using, for example Figure 3 The WOS method is shown for simulation. In each jump step, with the starting point x as the center of the sphere, the Green's function of the sphere indicates that the position of the next point in the reflected Brownian motion is uniformly distributed on the surface of the sphere. Therefore, a series of walks or jumps on the sphere can be used to simulate the path of the reflected Brownian motion. To speed up the WOS process, the radius of the sphere for each jump is set to the maximum value, i.e., the shortest distance from the current point to the boundary. When the point of the reflected Brownian motion is near the Dirichlet boundary, the point will be absorbed by the boundary, and the reflected Brownian motion will terminate; when the point of the reflected Brownian motion is near the Neumann boundary or Robin boundary, the local time needs to be calculated, and then the point will be reflected inward by the boundary.

[0070] In the heat source region, the WOG method is needed to realize reflected Brownian motion. This requires uniformly discretizing the heat source region into rectangular cuboid mesh elements, such as... Figure 4a As shown, the mesh needs to be fine enough to reasonably assume that points within each mesh cell have the same temperature and power density. Particles jump between the central nodes of these cuboid mesh cells, as... Figure 4b As shown, the jump probability in each direction is determined by the thermal conductance of adjacent grid nodes. The probability of jumping from grid node i to grid node j is...

[0071]

[0072] Where g(i,j) is the thermal conductance between nodes i and j, and d(i) is the total number of grid nodes adjacent to node i. Completing this jump will also contribute to the temperature, participating in the path integral of (6):

[0073]

[0074] Where p(i) is the power generated by the grid with center node i.

[0075] Therefore, substituting (10) into (6) and discretizing, the probabilistic solution for temperature can be modified as follows:

[0076]

[0077] Where N represents the total number of random walk paths, NP represents the maximum number of steps in a path, j' represents each step in the path, p represents the number of steps in the path within the heat source region, and j represents the number of steps in the path that collide with Neumann or Robin boundaries. The sum of the discrete forms of (4) and (5) represents the total temperature contribution under the three types of boundary conditions, dL(t). j ) is the localtime increment of the boundary corresponding to the collision point.

[0078] For those with N g The heat source region is a cubic grid, and since particles only jump to the center point of these grids, the power p(X) tp ) and thermal conductivity g t (X tp The grid number k = 1...N is only a function of the location of the heat source region. g to replace (11) can be rewritten as

[0079]

[0080] Since the power p(k) of the heat source is the only variable term in thermal management, (14) can be rewritten in vector form.

[0081]

[0082] If k' hotspots H1, H2, ..., H' are needed simultaneously k' The temperature, the matrix form of equation (15) is:

[0083]

[0084] Taking into account the fixed geometry and boundary conditions of the chip, the processing of matrix A and vector A can be completed offline. The pre-calculation is then performed, and then, during the real-time online operation of chip thermal management, the temperature prediction can be quickly completed through simple matrix multiplication and addition operations of equation (18).

[0085] It was observed that the pPIRW predictor (18) has a very similar structure to the predictor based on the compact thermal model (CTM), such as Figure 5As shown, the equivalence of the two can be proven through numerical experiments. In the construction of the CTM predictor, the thermal resistance matrix G needs to be inverted. When the chip is large and complex, this inversion process consumes a significant amount of time and space resources. The A matrix in the pPIRW predictor (18) based on the stochastic method is equivalent to G. -1 The columns corresponding to the hot spots are shown in the image. Figure 5 As shown, it is equivalent to solving a partial inverse of the large-scale thermal resistance matrix G.

[0086] In step 3 of this invention, a line from hotspot H is initialized. i The starting random walk path, i.e., the vector 'a' initialized in the pre-computation process. and boundary contribution value Record array.

[0087] For each hotspot H i Execute N random walking paths starting from the hotspot H, and in each path, move the hotspot H... i As the initial point X of the random algorithm t＝0 Initialize the pre-computation process. For the i-th hotspot H i ,make

[0088] a i =0, (20)

[0089] For from H i Let the j-th path be the starting point.

[0090]

[0091] They will be continuously updated in subsequent steps, ultimately based on (19) for matrix A and vector. The construction of.

[0092] In step 4 of this invention, based on the current point X t At the given location, perform the corresponding random walk procedure.

[0093] Step 4.1 If the current point X t In the passive region, a WOS random walk is performed

[10] .

[0094] Step 4.2 If the current point X t When the distance to the first type of Dirichlet boundary is less than ò1, then according to (13), the boundary value φ1 at the first type of boundary point closest to the current point is used to update the boundary. And terminate the random walk directly.

[0095] Step 4.3 If the current point X tD near the second type Neumann boundary ò Within the region, the boundary is calculated locally and updated according to (13). Specifically, in order to calculate local time

[20] , a strip region D with a width of ò is constructed near the Neumann boundary. ò ,like Figure 6 As shown, in D ò Within the region, the particle follows a fixed-radius WOS scheme, r = Δx or 2Δx. Assume the particle first enters D at x1. ò It then performs a WOS with a fixed radius Δx. When the particle jumps to a position less than Δx from the Neumann boundary (e.g., x2), the WOS radius is adjusted to 2Δx. The particle may then jump out of region Ω to x3, at which point it needs to be pulled back to the nearest position on the boundary, x4, and the WOS process continues from x4. According to (13), in the update... The boundary value φ2(x4) at position x4 on the boundary needs to be used.

[0096] For a certain boundary collision point Its corresponding boundary local time increment is

[11]

[0097]

[0098] in, X represents tj The corresponding boundary local time index increment is defined as

[0099]

[0100] Where, n j,out -n j,in Indicates local time from entry to exit (D) ò The coefficient increment of the entire process in the region

[11] , m j This represents the number of boundary collisions during this process.

[0101] Step 4.4 If the current point X t D near the third type Robin boundary ò Within the region, calculate the local time and update it. and local time and The update method is the same as in step 4.3, but in addition, the Feynman-Kac functional also needs to be updated. Update. The discrete form of equation (7) is:

[0102]

[0103] Among them, due to It is a negative value, and it only depends on the Robin boundary, therefore The value is updated only on each collision with a Robin boundary, decaying from 1 to 0. It's worth noting that... It will affect all path integral terms in (4)-(6).

[0104] Step 4.5 If the current point X t In the active region, a WOG random walk is performed. Based on the grid position number of the current point, if the current point is the k-th grid node in the heat source region, then a(k) is updated as follows:

[0105]

[0106] Then perform a WOG random walk.

[0107] In this invention, Figure 7 Show a pre-computed specific path that starts from x0 and... The path collides with the Robin boundary, passes through nodes 3 and 6 in the heat source region, and is eventually absorbed by the Dirichlet boundary. When passing through nodes 3 and 6, a(3) and a(6) are updated respectively. The physical meaning of this path is that the heat generated by nodes 3 and 6 is absorbed by the Robin boundary points. Therefore, a(k) is the probability that the heat generated by the k-th grid node will diffuse to the Robin boundary. In fact, the path of the random walk is entirely determined by the geometry, or more precisely, by the Green's function of the domain, and can therefore be calculated and stored offline in advance.

[0108] In step 5 of this invention, it is determined that X will jump after the jump. t+1 The position of X t+1 If the random walk is located on the Dirichlet boundary, the random walk ends, and proceed to step 6. If the random walk is not absorbed by the boundary, check the current... The size, if The random walk ends, and proceeds to step 6. If the number of random steps reaches the maximum number of jumps NP, the random walk ends, and proceeds to step 6. Otherwise, the jump to X is... t+1 As the current point X t Proceed to step 4 to continue the random walk.

[0109] In step 6 of this invention, the current random walk path terminates, and it is determined whether all the temperature hotspots to be measured have completed a random walk along N paths. If so, proceed to step 7; otherwise, proceed to step 3 to continue the random walk until all a... i and Pre-calculation of i = 1, 2, ..., k', j = 1, 2, ... N.

[0110] In step 7 of this invention, the sum of vector a obtained after executing N paths for each hotspot is calculated. The average value of all hotspots is calculated, and the pre-computed values ​​of all hotspots are assembled into matrix A and vector A.

[0111] For each hotspot, calculate the sum of vector a obtained after executing N paths. The average value, that is, for the i-th hotspot, is:

[0112]

[0113] According to (19) the matrix A and vector After the pre-calculation, the temperature prediction can be completed according to formula (18).

[0114] In step 8 of this invention, the pre-calculated matrix A and vector... This method is applied to active thermal management for real-time temperature prediction, and employs task mapping and DVFS measures to control chip hotspot temperatures to remain above a threshold. The specific process is as follows: Figure 8 As shown, the offline part consists of the matrix A and vector related to the hotspots obtained in steps 1-7. The pre-computation process, in its online component, is divided into two stages: task mapping and DVFS. In chip thermal management, multiple cores can function as a cluster, sharing the same voltage and frequency, and can be assigned tasks and perform DVFS as a whole.

[0115] During the task mapping phase, tasks equal to the number of idle clusters are retrieved from the task queue and sorted in descending order of average power. They are then assigned to the idle clusters furthest from the hottest clusters and with the lowest average temperature. During the DVFS tuning phase, the voltage and frequency levels of each cluster are adjusted progressively from highest to lowest until there are no thermal conflicts or all clusters are running at the lowest voltage and frequency levels. In each tuning, after determining the power of each computing core at the corresponding voltage and frequency levels, the online part (18) of the pPIRW predictor is invoked to obtain the hotspot temperature of all computing cores, calculate the number of thermal conflicts for each cluster, and reduce the voltage and frequency levels of the cluster with the most thermal conflicts to minimize the number of thermal conflicts.

[0116] The advantages of the method of the present invention are:

[0117] (1) By separating the variable power from other invariant components in the PIRW solver, the pPIRW method is derived, which yields an analytical solution for the steady-state temperature. As an accurate, fast, and physically interpretable temperature predictor, pPIRW is suitable for active thermal management of chips. The attenuation in the heat source region is pre-calculated offline. With fixed terms, temperature prediction can be obtained online via microsecond-level vector multiplication during real-time chip operation. This method provides probabilistic solutions for local regions, thus allowing direct calculation of the temperature value at a point within the domain without the need for global discretization of the region.

[0118] (2) After careful comparison of the temperature expressions of pPIRW and CTM, it was found that pPIRW is essentially a partial inverse of the large-scale thermal resistance matrix G in CTM calculated by a random method. This relationship has been numerically verified.

[0119] (3) pPIRW is a non-discretization method with constant time complexity. It does not require constructing a large-scale linear equation system or calculating matrix inverses, making it particularly suitable for engineers to handle large-scale or multi-scale problems with complex geometries. This method belongs to the stochastic method and has high parallelism. The stochastic paths are independent of each other and are easily parallelized, making this method very suitable for large-scale parallel computing.

[0120] This invention addresses the proactive thermal management problem during chip operation, proposing a pre-computed path integral random walk method as a powerful tool for rapid online temperature prediction and management. In the offline stage, it addresses large-scale thermal simulation problems involving complex structures, specifically solving mixed boundary value problems of the Laplace and Poisson equations. The method offers a significant speed advantage when applying massively parallel computing techniques. Furthermore, this method can be extended to other engineering fields for rapid field analysis and prediction problems involving Laplace and Poisson equations with mixed boundary conditions. Attached Figure Description

[0121] Figure 1 The algorithm flowchart of this method;

[0122] Figure 2 Schematic diagram of the definition of steady-state thermal analysis problem;

[0123] Figure 3 A diagram illustrating the WOS method;

[0124] Figure 4a Schematic diagram of heat source region discretization using the WOG method;

[0125] Figure 4b WOG method jump diagram;

[0126] Figure 5A schematic diagram illustrating the relationship between pPIRW and the CTM predictor;

[0127] Figure 6 The path is near the Neumann and Robin border, D. ò A walking map within the area;

[0128] Figure 7 pPIRW path diagram;

[0129] Figure 8 Active thermal management flowchart;

[0130] Figure 9 A schematic diagram of the chip's thermal model and the applied power.

[0131] Figure 10 Implementation example 2: Parallelism diagram. Detailed Implementation

[0132] The method of the present invention will now be described through specific examples of its implementation.

[0133] Experimental conditions settings:

[0134] like Figure 9 As shown, the test case structure consists of a three-layer thermal conductivity stack. The gray area represents the heat source region that continuously generates heat, simulating the die in the packaged chip. The top and bottom cuboid regions are passive regions, simulating the heat sink, dielectric layer, PCB, etc., in the packaged chip. The top surface represents the main heat transfer path through the heat sink, modeled with a convection coefficient h = 8700 W / m. 2 ·℃, ambient temperature T amb The Robin boundary is at 20°C. The bottom surface represents a secondary heat transfer path through the PCB, which can be modeled as isothermal (Dirichlet, T = 70°C), adiabatic (Neumann), and convective (Robin, h = 8700 W / m²). 2 ·℃,T amb With boundary conditions of 20℃, there is virtually no heat conduction on the four sides, so it is modeled as an adiabatic boundary. The power diagram applied to the heat source region is shown in the top view. Figure 9 As shown, 4core-1 and 4core-2 are power graphs of a randomly generated 4-core chip, while 16core is a power graph generated by the POWER6 microprocessor. 16core is a power graph generated through real-time simulation for thermal management of a 16-core chip.

[0135] The three-layer stacked region is 2cm in length and width, and its thickness and thermal conductivity are shown in the table below:

[0136] Table 1: Thickness and thermal conductivity of the three-layer structure

[0137] Thickness / um Thermal conductivity (W / K·m) Top floor 500 395 Mid-level heat source region 100 125 bottom layer 500 395

[0138] In this method, the number of paths N = 5e3, and is reduced to 1.5e3 in cases without Dirichlet boundaries to speed up the process. The maximum number of steps NP = 6e6, the path termination error tolerance ε = 5e-4 meters, and the distance near the boundary D ò The internally fixed WOS radius Δx = 5e-7 meters, ò = 1.36Δx, and the absorbing boundary band width ò1 = 1e-8 meters.

[0139] Implementation Example 1

[0140] The objective of this example is to compare the accuracy and speed of the solution results obtained by the patented method pPIRW with those of commercial software ANSYS, academic open-source tools 3D-ICE and HotSpot, and the stochastic method PIRW.

[0141] The top surface of the Mix4 and Stack-1 tests has a Robin boundary condition, and the bottom surface has a Dirichlet boundary condition. Power maps of 4core-1 and 4core-1+4core-2 stacks are applied respectively in the heat source region. The top and bottom surfaces of the R4 and Stack-2 tests both have Robin boundary conditions, and power maps of 4core-1 and 4core-1+4core-2 stacks are applied respectively in the heat source region.

[0142] Table 2 shows the speed and accuracy comparison results of the pPIRW method with ANSYS, 3D-ICE, HotSpot, and PIRW. In all four examples, the errors of pPIRW and PIRW methods are within 0.152℃. In examples containing Dirichlet boundaries (Mix4 and Stack-1), PIRW achieves a speedup of 136-209 times and a memory reduction of 8.53-11.1 times compared to 3D-ICE; in examples without Dirichlet boundaries, this speed advantage degrades to 1-1.53 ​​times. pPIRW only requires 3.46-6.92us to complete temperature prediction, achieving a speedup of 2.5×10⁻⁶ times compared to PIRW. 4 -6.7×10 6 A speedup of 1x and a memory reduction of 36.6-42.5x.

[0143] Table 2: Comparison of speed and accuracy of the patented method pPIRW with ANSYS, 3D-ICE, HotSpot, and PIRW

[0144]

[0145] *The values ​​in parentheses represent the time and memory calculatedoffline by pPIRW.

[0146] Implementation Example 2

[0147] The goal of this example is to test the parallel efficiency of the method of this invention. The settings of the Mix4 example are the same as those of Example 1. Figure 10 The experimental results shown indicate that the randomization method of this application can maintain an approximately linear speedup when there are 80 CPU cores (maximum number of cores per machine).

[0148] Implementation Example 3

[0149] The goal of this example is to numerically verify the consistency between the pPIRW method and the CTM method of this patent. The thermal resistance matrix is ​​constructed using MatEx

[16] and inverted using the Eigen library

[22] . Figure 9 In the thermal model, a 1.6mm × 1.6mm × 1.1mm block is vertically cut. It still has a three-layer cubic stacked structure. The top and bottom layers are non-heat source areas with a thickness of 0.5mm, and the middle layer is a heat source area with a thickness of 0.5mm. The heat source area is divided into a grid of 16 × 16 × 5 = 1280. According to Figure 5 Each point in pPIRW will have a 1280-dimensional a' vector, which is compared with the corresponding point in the CTM method from the inverse of the thermal resistance matrix G. -1 The similarity of vector a extracted from it.

[0150] As shown in Table 3, the relative error between vector a' and vector a is only 2.34%. Temperature prediction requires vector multiplication, so its error depends on the error of the sum of elements of vectors a' and a, which is only 0.18%. Experimental results verify the consistency between the method of this patent and CTM. Therefore, pPIRW is essentially solving for a partial inverse of the large-scale thermal resistance matrix in CTM.

[0151] Table 3 shows the error between the pre-calculated vector a' and the partial inverse a of the thermal resistance matrix in the CTM using the method of this patent.

[0152]

[0153] Implementation Example 4

[0154] The goal of this example is to compare the effectiveness of the pPIRW method in active thermal management.

[0155] To simulate a more realistic chip package, the thickness, width, and thermal conductivity of the bottom heat source region were modified to 100µm, 1.6cm×1.6cm, and 50W / K·m, respectively. The ambient temperature was set to 45℃, and the top surface was a Robin boundary (h = 12500 or 10000W / m). 2 ·℃), using Figure 9 The power plot of the 16 cores is shown. The assigned task is from the PARSEC v2.1 test set, using four V / f levels. The transient simulation results from ANSYS are used as the reference solution.

[0156] The effects of active thermal management methods using different predictors are shown in Table 4. The average error of the pPIRW and CTM methods in predicting temperature is less than 0.27℃, and the maximum error is less than 1.15℃; while the average error of the RM method in predicting temperature is greater than 3.82℃. The pPIRW and CTM methods can basically eliminate all thermal conflicts and keep the chip operating at the highest possible frequency; moreover, they are both based on physical models and do not require access to thermal sensors, controlling single-step thermal management at the microsecond level.

[0157] As shown in Table 4, in the offline pre-computation stage, this method achieves a speedup of 29.2 times and a memory saving of 634 times compared to CTM. The speed and memory improvement stem from the fact that the patented method pPIRW avoids the global discretization, linear equation system construction, and time-consuming large matrix inversion required by deterministic methods (such as CTM). CTM requires constructing a large-scale thermal resistance matrix and performing matrix inversion, with a time complexity of O(n^2). 3 ), where n is the degree of freedom. Since the inverse of a sparse matrix is ​​usually not sparse, directly inverting and storing a large sparse thermal resistance matrix G using CTM in practical large-scale thermal analysis is very time-consuming and space-intensive. Especially when only the temperature of a portion of the chip's hot spots is needed, only G- is actually required. 1 For some elements of G, performing a complete inversion of G becomes extremely tedious and inefficient. However, the patented method pPIRW can solve for the temperature of only some hotspots, with a time complexity of O(cm), where m is the time to calculate the temperature of a single hotspot using a random method, and c is the number of hotspots. Therefore, this patented method avoids time-consuming global discretization, linear system construction, and matrix inversion. In large-scale chip and complex geometry problems, the advantages of pPIRW in terms of speed, memory usage, lack of discretization, and parallelism will be even more pronounced.

[0158] Table 4: Effects of Active Thermal Management Methods Using Different Predictors

[0159]

[0160] *lnclude he time to read theral sensors trom

[34] 。

Claims

1. A chip thermal management method based on a pre-computed random method steady-state thermal analysis solver, used to realize a high-precision, high-speed temperature predictor in chip active thermal management; characterized in that, The thermal management problem of a chip is divided into an offline thermal analysis problem and an online dynamic thermal management optimization problem; In the offline thermal analysis stage, the chip thermal analysis problem is transformed into a Poisson equation solving problem with mixed boundary conditions. Based on the Feynman-Kac formula, the WOS method of standard reflected Brownian motion is used to absorb or reflect the first, second and third types of boundaries respectively, and the probability solution of the hotspot temperature in the chip steady-state thermal analysis problem is obtained. The power term is decoupled from the random walk path solution, and the changing power term in the chip thermal management problem is extracted from other unchanged parts through pre-computation by random walk, and the pre-computed matrix A and vector and the temperature prediction formula in the form of fast matrix multiplication and addition is obtained; In the online dynamic thermal management stage, the pre-computed matrix A and vector The application is applied to the active thermal management method to realize real-time prediction of temperature, and uses task mapping and dynamic voltage frequency adjustment (DVFS) measures to control the chip hotspot temperature not to exceed the set temperature threshold; specifically includes: Input parameters: region Ω and its boundaries Geometric information, material thermal conductivity parameters, ambient temperature T amb ,boundary The three types of boundaries Γ i Boundary conditions φ i Let i = 1, 2, 3, N be the number of Monte Carlo sampling paths, NP be the maximum number of path jumps, ε be the path termination error tolerance, and k' hotspot locations H1, H2, ..., H k' D near the boundaries of the second and third categories ò The thickness of the region ò and the jump radius Δx of the random walk inside, the width ò1 of the absorption band near the first-class boundary, and the power diagram of the heat source region over time. Output parameters: the steady-state temperature of the chip hot spot is predicted, and the task mapping and the dynamic voltage and frequency scaling (DVFS) level are adjusted according to the predicted temperature; The method comprises the following steps: Step 1: the thermal management problem of a chip is divided into a thermal analysis problem and an online dynamic thermal management optimization problem, and the chip thermal analysis problem is further converted into a Poisson equation solving problem with mixed boundary conditions; Step 2: Based on the Feynman-Kac formula, the probability solution of hotspot temperature of the chip steady-state thermal analysis problem is obtained, the power term is decoupled from the random walk path solution, and the temperature prediction formula in the form of matrix multiplication and addition at k' hotspots (H1, H2,..., H k' ) is obtained; Step 3: Initialize a random walk path starting from hotspot Hi, i.e. initialize the a vector of the precomputation process, and the boundary contribution value the record array; Step 4: According to the position of the current point Xt, execute the corresponding random walk process, use WOS random walk in the passive region; update near the first type of Dirichlet boundary and directly terminate the random walk; calculate the local time and update near the second type of Neumann boundary calculate the local time and update near the third type of Robin boundary and In the active region, execute the WOG random walk, and update the corresponding position elements of the a vector according to the grid position number where the current point is located; Step 5: Determine Xt after the jump +1 The position, if Xt +1 Located on the absorption boundary or If the maximum number of jumps in the path is reached (NP), the random walk ends and proceeds to step 6; otherwise, the jump to the next step (Xt) is completed. +1 As the current point Xt, proceed to step 4; Step 6: it is judged whether the k' hot spots have completed N random walking paths, if yes, the step 7 is proceeded, otherwise, the step 3 is proceeded to continue random walking; Step 7: Calculate the average of the a-vectors obtained for each hotspot after performing N paths and assemble the pre-computed values of all hotspots into matrix A and vector b ​ Step 8: Pre-computed matrix A and vector The real-time prediction of temperature is applied to the active thermal management method, and the task mapping and DVFS measures are used to control the chip hotspot temperature not to exceed the threshold value.

2. The method of claim 1, wherein In step 1, the thermal management problem of a chip is divided into a thermal analysis problem and an online dynamic thermal management optimization problem, and the chip thermal analysis problem is further converted into a Poisson equation solving problem with mixed boundary conditions; The main heat dissipation path of the chip is through a heat sink, so the top surface of the chip is modeled as a Robin boundary; the secondary heat transfer path is through a PCB, so the bottom surface of the chip can be modeled as a Dirichlet, Neumann or Robin boundary according to the actual situation; the four side surfaces are usually modeled as adiabatic Neumann boundaries because there is almost no heat conduction; In the mathematical model of three-dimensional steady-state heat conduction problem, the three-dimensional region of steady-state heat analysis is Ω, the temperature field T(x) in the region Ω satisfies Poisson equation, the boundary of the region Ω is a mixed boundary, that is, contains Dirichlet boundary Γ1, Neumann boundary Γ2 and Robin boundary Γ3, or only one or two of them; given three types of boundary conditions, then the field differential equation and boundary conditions are: where T(x) is the temperature in the domain Ω, p d (x) is the heat generation rate, k(x) is the thermal conductivity of the material; on the Dirichlet boundary Γ1, the temperature T(x) is known as φ1(x); on the Neumann boundary Γ2, the outward normal derivative of the temperature T at the boundary is known as φ2(x); in steady-state heat analysis, it is usually an adiabatic boundary condition, i.e. Robin boundary, also known as convective boundary, corresponds to the interface that exchanges heat with the surrounding fluid, such as the surface of cooling fans and heat sinks; the parameter c(x) in the Robin boundary condition is: Robin boundary, also known as convective boundary, corresponds to the interface that exchanges heat with the surrounding fluid, such as the surface of cooling fans and heat sinks; the parameter c(x) in the Robin boundary condition is: where h(x) is the heat convection coefficient at Robin boundary, which is usually where T amb is the known ambient temperature; The chip active thermal management process is expressed as an optimization problem: where the optimization objective is to maximize speed, i.e. maximize the average voltage or frequency of the computing cores, while keeping the hotspot temperature T at a threshold value T threshold The following.

3. The method of claim 1 wherein, In step 2, the probability solution of hotspot temperature of the chip steady-state thermal analysis problem is obtained based on the Feynman-Kac formula, the power term is decoupled from the path solution of random walk, and a temperature prediction formula in the form of matrix multiplication and addition at k' hotspots (H1, H2,..., H k' ) is obtained. Considering the Poisson equation (1) with mixed boundary conditions, according to the Feynman-Kac formula, the probability solution of the temperature value of a point x in the region Ω is: where E x is the mathematical expectation, Xtis the reflected Brownian motion, is the time when the reflected Brownian motion reaches the Dirichlet boundary, is the position corresponding to the time; is the Feynman-Kac functional, L(t) is the boundary local time corresponding to the reflected Brownian motion near the Neumann boundary, Robin boundary, where D ò is the strip region near the boundary with width ò, is the indicator function in the region D ò ; the integral terms in (4)-(6) represent the contributions of the right-hand side of the Poisson equation to the temperature at the point under test on the Dirichlet boundary, Neumann and Robin boundary, respectively. The reflected Brownian motion is usually simulated by using the WOS method; in each jump step, the starting point x is taken as the center of a sphere, and according to the Green function of the sphere, the position of the next point of the reflected Brownian motion is uniformly distributed on the surface of the sphere; therefore, a series of walks or jumps on the sphere can be used to simulate the path of the reflected Brownian motion; in order to speed up the process of WOS, the radius of the sphere in each jump is set to the maximum value, that is, the shortest distance from the current point to the boundary; when the point of the reflected Brownian motion is near the Dirichlet boundary, the point is absorbed by the boundary, and the reflected Brownian motion is terminated; when the point of the reflected Brownian motion is near the Neumann boundary or the Robin boundary, the local time needs to be calculated, and then the point is reflected to the inside by the boundary; In the heat source region, the WOG method is used to realize the reflected Brownian motion; it is necessary to uniformly discretize the heat source region into rectangular cuboid grid cells, and the grid needs to be fine enough to reasonably assume that the points in each grid cell have the same temperature and power density; particles jump between the center nodes of these rectangular cuboid grid cells, and the jump probability in each direction is determined by the thermal conductivity of the adjacent grid nodes; the probability of jumping from grid node i to grid node j is Where g(i,j) is the thermal conductivity between nodes i and j, and d(i) is the total number of grid nodes adjacent to node i; completing this jump also contributes to the temperature, and the path integral participating in (6) is Where p(i) is the power generated by the grid with the center node i. Therefore, plugging (10) into (6) and discretizing, the probabilistic solution of temperature can be modified as where N denotes the total number of random walks, NPthe maximum number of steps of the path, j' represents each step in the path, p the number of steps in which the path is inside the heat source region, and j the number of steps in which the path hits the Neumann or Robin boundary; is the sum of (4) and (5) in discrete form, representing the total sum of temperature contributions under the three types of boundary conditions, dL(t j ) is the boundary local time increment corresponding to the collision point; For a heat source region with N g cubic grids, since the particles will only hop to the center points of these grids, the power and heat conductance are only functions of the heat source region location, so we can replace g (11) with​ Since in thermal management, the power p(k) of heat source is the only changing item, (14) can be rewritten in vector form as If k' hotspots H1, H2, ..., H' are needed simultaneously k' The temperature, the matrix form of equation (15) is: Where the matrix A and vector b can be pre-computed in the off-line part considering the geometry and boundary conditions of the chip fixation, and then the temperature prediction can be quickly completed by the simple matrix multiplication operation of equation (18) in the real-time on-line operation of chip thermal management. Where the matrix A and vector b can be pre-computed in the off-line part considering the geometry and boundary conditions of the chip fixation, and then the temperature prediction can be quickly completed by the simple matrix multiplication operation of equation (18) in the real-time on-line operation of chip thermal management. The pPIRW predictor (18) and the compact thermal model based predictor have similar structures, and their equivalence can be proved by numerical experiments. In the construction of the CTM predictor, the inverse of the thermal resistance matrix G is needed, which consumes a large amount of time and space resources when the chip scale is very large and the structure is complex. The A matrix in the pPIRW predictor (18) based on the random method is equivalent to the part of the columns of the part of the rows corresponding to the hot spots in the G -1 matrix, so it is equivalent to solving the partial inverse of the large-scale thermal resistance matrix G.

4. The method of claim 1 wherein, In step 3, a random walk path is initialized from the hotspot Hi, i.e. the a vector of the precomputation process is initialized, and the boundary contribution value the record array; For each hotspot Hi, N random walk paths are performed from the hotspot, and in each path, the hotspot Hi is taken as the initial point of the random algorithm Xt=0, and the pre-computation process is initialized. For the i-th hotspot Hi, let a i ＝0, (20) For the j-th path from Hi, let They will be updated in the subsequent steps and eventually used for the construction of the matrix A and the vector b according to (19) .

5. The method of claim 1 wherein, In step 4, according to the position of the current point Xt, the corresponding random walk process is performed; Step 4.1 If the current point Xt is in the source-free area, WOS random walk is performed; Step 4.2 If the current point Xtis near the Dirichlet boundary of the first kind, i.e., the distance to the Dirichlet boundary is less than 0i, then update the boundary value f at the current point according to (13) with the boundary value f at the nearest boundary point of the first kind, and directly terminate the random walk. and directly terminate the random walk. Step 4.3 If the current point Xt is near the second kind of Neumann boundary D ò Within the region, the boundary is calculated locally and updated according to (13). Specifically, to calculate local time, a strip-shaped region D with a width of ò is constructed near the Neumann boundary. ò In D ò Within the region, the particle follows a fixed-radius WOS scheme, r = Δx or 2Δx; assuming the particle first enters D at x1. ò It then performs WOS with a fixed radius Δx; when the particle jumps to a position less than Δx from the Neumann boundary, such as x2, the WOS radius is adjusted to 2Δx; next, the particle may jump out of region Ω to x3, at which point it needs to be pulled back to the nearest position on the boundary x4, and the WOS process continues from x4; according to (13), in the update At this time, the boundary value φ2(x4) at position x4 on the boundary needs to be used; For a certain boundary collision point The corresponding boundary local time increment is wherein denotes the corresponding boundary local time indicator increment, defined as where n j,out - n j,in represents the coefficient increment of the local time from entering to leaving the D ò region, m j is the number of boundary collisions in this process; Step 4.4 If the current point Xtis in the region D near the third type of Robin boundary, then compute the local time and update ò the local time and the local time and the local time and in the same way as in Step 4.3, but in addition, the Feynman-Kac functional needs to be updated; the discrete form of equation (7) is: whereas is a negative value, which is only related to Robin boundary, thus is only updated in the path at each collision with Robin boundary, and its value decays from 1 to 0 gradually; it is worth noting that will have an impact on all path integral terms in (4)-(6); Step 4.5 If the current point Xt is in the source area, WOG random walk is performed, according to the grid position number where the current point is located, if the current point is the k-th grid node in the heat source area, update a(k) as: Then WOG random walk is performed.

6. The method of claim 1 wherein, The average of the a vectors obtained after performing N paths for each hotspot in step 7 is calculated, and the pre-computed values of all hotspots are assembled into matrix A and vector b. ​ For each hotspot, compute the average of the a-vectors it gets after executing N paths, i.e., for the i-th hotspot, have ​ After pre-computation of the matrix A and the vector according to (19), the temperature prediction can be done according to formula (18).

7. The method of claim 1 wherein, The pre-computed matrix A and vector The real-time prediction of temperature is applied to the active thermal management method, and the task mapping and DVFS measures are used to control the chip hotspot temperature not to exceed the threshold value; the offline part is the pre-computation of the matrix A and vector of the hotspot completed by steps 1-7, and the online part is divided into two stages of task mapping and DVFS; in chip thermal management, multiple cores can be a cluster, which share the same voltage and frequency, and can be allocated tasks and DVFS as a whole; In the task mapping stage, the same number of tasks as the number of idle clusters are taken out from the task queue, and they are sorted in descending order of average power, and are allocated to the idle clusters far from the hottest cluster and with the lowest average temperature in turn; In the DVFS adjustment stage, the voltage and frequency levels of each cluster need to be adjusted from the highest to the lowest step by step until there is no thermal conflict or all clusters run at the lowest voltage and frequency level; In each adjustment, after determining the power of each computing core of the task at the corresponding voltage and frequency level, the online part (18) of the pPIRW predictor is called to obtain the hotspot temperature of all computing cores, the number of thermal conflicts of each cluster is calculated, and the voltage and frequency level of the cluster with the most thermal conflicts are reduced to minimize the number of thermal conflicts.

Images