A high-precision random method three-dimensional chip steady-state thermal analysis acceleration method based on GPU parallelism and fitting

Abstract

The present application belongs to the field of integrated circuit computer aided design / electronic design automation (Computer Aided Design / Electronic Design Automatic CAD / EDA), and particularly relates to solving the problem of steady-state thermal analysis under mixed boundary conditions in three-dimensional chips. The method uses GPU parallel and fitting technology to accelerate the solving of mixed boundary value problems (mixed BVPs) of Poisson equation in the steady-state thermal analysis of three-dimensional chips based on the path integral random walk (PIRW) method, i.e. the partial differential equations (PDE). The method can fully utilize the advantages of the large amount of GPU computing resources, effectively accelerate the processing of problems containing the second and third boundary conditions, and particularly for the slow decay problem of the Feynman-Kac functional related to the third Robin boundary, the parallel efficiency is further improved by using the fitting technology.

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 the solution of steady-state thermal analysis problems under mixed boundary conditions in three-dimensional chips. This method utilizes GPU parallelism and fitting techniques to effectively accelerate the processing of the problem of slow decay of the Feynman-Kac functional related to Robin boundaries, further improving the parallel efficiency. Background Technology

[0002] In existing technologies, with the continuous shrinking of feature sizes and the development of integrated circuit design, Chiplet / 3D ICs have become an important technological direction. However, thermal issues significantly impact the long-term reliability and performance of integrated circuits. Therefore, accurate and efficient thermal analysis is indispensable to ensure that chip design and operation meet temperature constraints. Meanwhile, GPUs, due to their abundant parallel computing resources, are widely used and have proven highly efficient in solving many electronic design automation (EDA) problems.

[0003] Existing thermal simulation methods can be divided into three categories: deterministic methods, Green's function methods, neural network-based methods, and stochastic methods.

[0004] Deterministic methods include the finite element method (FEM)[1],[2] and the finite difference method (FDM)[3],[4]. The finite element method discretizes a continuous region into interconnected elements and solves partial differential equations using variational principles. The finite difference method is based on the differential approximation at the grid points and constructs a system of algebraic equations by replacing the derivative with the difference quotient. Although deterministic methods are highly accurate and adaptable, solving large-scale problems is very time-consuming because it generally involves solving a large number of linear equation systems.

[0005] The Green's function method [5] and [6] transforms the calculation of the field into the convolution of the Green's function and the source, but the Green's function usually does not have an analytical form in practical problems, and the numerical solution needs to be obtained through deterministic methods.

[0006] Neural network-based methods directly fit accurate labels [7] or use physical information loss functions [8] to predict the estimated temperature value of any point through the network. However, they are difficult to handle geometric changes, and the network training time is long and the error range cannot be estimated.

[0007] As a potential solution for high-precision hotspot analysis

[12] ,

[13] ,

[15] , the stochastic method has gradually shown its unique advantages in the field of chip thermal simulation. Based on the Feynman-Kac formula [9], the stochastic method obtains the single-point solution of the partial differential equation by releasing many particles from the point to be solved to perform Brownian motion and calculating the expected value of their path integral. The stochastic method has been widely used for capacitance extraction

[10] and has been widely implemented in commercial software. Path integral random walk (PIRW)

[13] is a high-precision solver based on the stochastic method. It uses rigorous and accurate local time calculation and Feynman-Kac functional.

[14] and

[16] handle mixed boundary conditions.

[0008] However, since the results obtained by the stochastic method generally have The variance of the path is given by the number of paths

[12] , which means that more than a thousand independent paths are needed to achieve high accuracy. If parallel technology can be used, the computational efficiency of the program can be greatly improved. However, since the third type of Robin boundary is common in chip thermal analysis problems, each random walk path usually needs to perform millions of jumps to ensure that the final temperature of all paths converges. This characteristic limits the parallel effect of the PIRW program. Even if all paths are parallelized, the final parallel execution time cannot be effectively reduced.

[0009] To improve the parallel performance of the PIRW program, this invention proposes a method to terminate random walks in advance before the temperature has fully converged, and to estimate the final convergence temperature by fitting. This method can further improve parallel efficiency and determine the timing for terminating random walks in advance.

[0010] This method addresses steady-state thermal analysis problems under mixed boundary conditions in 3D chips. The proposed GPU parallel acceleration method is a powerful tool for solving mixed boundary value problems of the Laplace and Poisson equations, offering a significant speed advantage for solving complex thermal problems in large-scale 3D chips. Furthermore, this method can be extended to field analysis problems of the Laplace and Poisson equations with mixed boundary conditions in other engineering fields.

[0011] References related to this invention;

[0012] [1]IncANSYS et al.Ansys.https: / / www.ansys.com / ,2018.

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

[0014] [3]Wei Huang,S.Ghosh,S.Velusamy,K.Sankaranarayanan,K.Skadron,andM.R.Stan.Hotspot:a compact thermal modeling methodology for early-stage vlsidesign.IEEE Transactions on Very Large Scale Integration(VLSI)Systems,14(5):501–513,2006.

[0015] [4]Federico Terraneo,Alberto Leva,William Fornaciari,Marina Zapater,and David Atienza.3d-ice 3.0:Efficient nonlinear mpsoc thermal simulationwith pluggable heat sink models.IEEE Transactions on Computer Aided Design ofIntegrated Circuits and Systems,41(4):1062–1075,2022.

[0016] [5]Hameedah Sultan and Smruti R.Sarangi.A fast leakage-aware green’sfunction-based thermal simulator for 3-d chips.IEEE Transactions on VeryLarge Scale Integration(VLSI)Systems,28(11):2342–2355,2020.

[0017] [6]Hameedah Sultan and Smruti R.Sarangi.Varsim:A fast and accuratevariability and leakage aware thermal simulator.In 202057th ACM / IEEE DesignAutomation Conference(DAC),pages 1–6,2020.

[0018] [7]Tianxiang Zhu,Qipan Wang,Yibo Lin,Runsheng Wang,and RuHuang.Fasttherm:Fast and stable full-chip transient thermal predictorconsidering nonlinear effects.In 2024 IEEE / ACM International Conference onComputer Aided Design(ICCAD),pages 1–9,2024.

[0019] [8]Liang Chen,Jincong Lu,Wentian Jin,and Sheldon X.-D.Tan.Fast fullchip parametric thermal analysis based on enhanced physics enforced neuralnetworks.In2023 IEEE / ACM International Conference on Computer Aided Design(ICCAD),pages 1–8,2023.

[0020] [9]EAN-PAULMORILLON.Numerical solutions of linear mixed boundaryvalue problems using stochastic representations.International journal fornumerical methods in engineering,40(3):387–405,1997.

[0021]

[10] Wenjian Yu,Hao Zhuang,Chao Zhang,Gang Hu,and Zhi Liu.Rwcap:A floating random walk solver for 3-d capacitance extraction of very-large scaleintegration interconnects.IEEE Transactions on Computer-Aided Design ofIntegrated Circuits and Systems,32(3):353–366,2013.

[0022]

[11] Kuangya Zhai,Wenjian Yu,and Hao Zhuang.2013.GPU-friendly floatingrandom walk algorithm for capacitance extraction of VLSI interconnects.InProceedings of the Conference on Design,Automation and Test in Europe(DATE’13).EDA Consortium,San Jose,CA,USA,1661–1666.

[0023]

[12] Yuan Liang,Wenjian Yu,and Haifeng Qian.A hybrid random walkalgorithm for 3-d thermal analysis of integrated circuits.In 201419th Asiaand South Pacific DesignAutomation Conference(ASP-DAC),pages 849–854,2014.

[0024]

[13] Longlong Yang,Cuiyang Ding,Changhao Yan,Dian Zhou,and Xuan Zeng.Ahigh-precision stochastic solver for steady state thermal analysis withfourier heat transfer robin boundary conditions.In 2022 IEEE / ACMInternational Conference On ComputerAided Design(ICCAD),pages 1–9,2022.

[0025]

[14] Cuiyang Ding,Yijing Zhou,Wei Cai,Xuan Zeng,Changhao Yan,A pathintegral Monte Carlo(PIMC)method based on Feynman Kac formula for electricalimpedance tomography,Journal of Computational Physics,Volume 476,2023,111862,ISSN 0021-9991,https: / / doi.org / 10.1016 / j.jcp.2022.111862.

[0026]

[15] Zhixuan Dong,Longlong Yang,Cuiyang Ding,Changhao Yan,Zhaori Bi,Sheng-Guo Wang,Dian Zhou,and Xuan Zeng.PPIRW:An efficient and accurate pre-calculation path integral random walk solver for steady state thermalsimulation with robin boundary conditions.IEEE Trans actions on ComputerAided Design of Integrated Circuits and Systems,pages 1–1,2024.

[0027]

[16] Vassilis G Papanicolaou.1990.The probabilistic solution of the third boundary value problem for second order elliptic equations.ProbabilityTheory and related fields 87,1(1990),27–77.

[0028]

[17] Pei Hsu.1984.Reflecting Brownian motion,boundary local time and the Neumann problem.Ph.D.Dissertation.Stanford University.

[0029]

[18] Yijing Zhou and Wei Cai.2016.Numerical solution of the Robinproblem of Laplace equations with a Feynman–Kac formula and reflecting Brownian motions.Journal ofScientific Computing 69,1(2016),107–121.

[0030]

[19] Yijing Zhou, Wei Cai, and Elton Hsu. 2017. Computation of the localtime of reflecting brownian motion and the probabilistic representation of the neumann problem. Communications in Mathematical Sciences 15, 1 (2017), 237–259. Summary of the Invention

[0031] The purpose of this invention is to propose a path integral stochastic method based on GPU parallelism and path fitting technology, based on the current state of existing technology, for GPU parallel acceleration to solve steady-state thermal analysis problems with mixed boundary conditions in three-dimensional chips.

[0032] In existing technologies, the results of randomized methods generally have... The variance of the path is given by the number of paths

[12] , which means that more than a thousand independent paths are needed to achieve high accuracy. If parallel technology can be used, the computational efficiency of the program can be greatly improved. However, since the third type of Robin boundary is common in chip thermal analysis problems, each random walk path usually needs to perform millions of jumps to ensure that the final temperature of all paths converges. This characteristic limits the parallel effect of the PIRW program. Even if all paths are parallelized, the final parallel execution time cannot be effectively reduced.

[0033] This method leverages GPU parallelism and fitting techniques to accelerate the solution of mixed boundary value problems (mixed BVPs) in 3D chip steady-state thermal analysis using the Path Integral Random Walk (PIRW) method. This method fully utilizes the abundant computing resources of GPUs to effectively accelerate the handling of problems with second and third type boundary conditions. Particularly for the challenge of slow decay of the Feynman-Kac functional related to the third type of Robin boundary, fitting techniques are used to further improve parallel efficiency.

[0034] Specifically, in this invention,

[0035] First, the steady-state thermal analysis problem of the 3D chip is transformed into a Poisson equation boundary value problem with mixed boundary conditions. Then, based on the Feynman-Kac formula, a probabilistic solution to the steady-state thermal analysis problem of the 3D chip is obtained using standard reflection Brownian motion. Third, all input parameters are transferred to the GPU, N GPU threads are started according to the required number of sampling paths N, and the program is initialized. Fourth, based on the current point X... t At the given location, a corresponding random walk process is executed in the GPU; random walks are performed using WOS and WOG in the passive and hot regions, respectively; for three different types of boundary conditions, direct termination, local time calculation, and local time calculation and update are adopted, respectively. Three different calculation strategies; Fifth, determine whether the number of jump steps in the current path exceeds an integer multiple of the jump step number δ, i.e., np≥kδ. If so, proceed to step six; otherwise, proceed to step four; Sixth, record the cumulative temperature T of the i-th path when the number of jump steps is kδ, calculated using the discrete form of the Feynman-Kac formula. i，kδ(x), if all N paths have completed at least kδ jumps, proceed to step 7; otherwise, proceed to step 4. Seventh, calculate the average total contribution of the third-type Robin boundary and heat source region to the cumulative temperature value when the number of jumps for the N paths is kδ. After adding this average value, the temperature is fitted within the temperature fitting window w. To determine whether the fitted temperature has converged, i.e., to determine... If yes, terminate the program and calculate the temperature T(x) at the point x to be determined as the output; otherwise, increment k by 1 and proceed to step four.

[0036] More specifically,

[0037] This invention proposes a path integral stochastic acceleration method based on GPU parallelism and temperature fitting for rapidly solving steady-state thermal analysis problems with mixed boundary conditions in 3D chips. The flowchart is shown below. Figure 1 As shown, it includes:

[0038] Input parameters: region Ω and its boundaries Geometric information, region boundaries The three types of boundaries Γ i and boundary conditions i = 1, 2, 3, near the first-kind absorbing boundary Thickness of the strip region Near the boundaries of Type II and Type III reflections Thickness of the strip region and Within the strip region, the WOS jump radius Δx, the location of the point x to be determined, the number of sampling paths N, the interval jump steps δ for initiating fitting, and the convergence criterion ε for judging the fitted value. th Temperature fitting window w;

[0039] Output: Temperature T(x) at point x to be determined.

[0040] The specific steps include:

[0041] Step 1: Transform the steady-state thermal analysis problem of the 3D chip into a boundary value problem of the Poisson equation with mixed boundary conditions;

[0042] Step 2: Based on the Feynman-Kac formula, and using standard reflection Brownian motion, obtain the probabilistic solution to the three-dimensional chip steady-state thermal analysis problem;

[0043] Step 3: Transfer all input parameters to the GPU, start N GPU threads according to the required number of sampling paths N, and initialize the program;

[0044] Step 4: Based on the current point X tAt the given location, a corresponding random walk process is executed in the GPU; random walks are performed using WOS and WOG in the passive and hot regions, respectively; for three different types of boundary conditions, direct termination, local time calculation, and local time calculation and update are adopted, respectively. Three different computational strategies;

[0045] Step 5: Determine whether the number of jump steps in the current path exceeds an integer multiple of the number of jump steps δ, i.e., np≥kδ. If yes, proceed to step 6; otherwise, proceed to step 4.

[0046] Step 6: Record the cumulative temperature T for the i-th path when the number of jump steps is kδ, calculated using the discrete form of the Feynman-Kac formula. i,kδ (x), if all N paths have completed at least kδ jumps, proceed to step 7; otherwise, proceed to step 4.

[0047] Step 7: Calculate the average total contribution of the third-type Robin boundary and heat source region to the cumulative temperature value when the number of jump steps for N paths is kδ. After adding this average value, the temperature is fitted within the temperature fitting window w. To determine whether the fitted temperature has converged, i.e., to determine... If yes, terminate the program and calculate the temperature T(x) at the point x to be determined as the output; otherwise, increment k by 1 and proceed to step 4.

[0048] In step 1 of this invention, the steady-state thermal analysis problem of a three-dimensional chip is transformed into a boundary value problem of the Poisson equation with mixed boundary conditions.

[0049] Specifically, the problem solved by the method of the present invention is as follows: Figure 2 As shown. The three-dimensional region for steady-state thermal analysis is... Within region Ω, the temperature field T(x) satisfies the Poisson equation, and the region boundary... The boundary conditions are mixed, meaning they include a first-kind Dirichlet boundary Γ1, a second-kind Neumann boundary Γ2, and a third-kind Robin boundary Γ3. Given these three types of boundary conditions, the field differential equations and boundary conditions are as follows:

[0050]

[0051] Where T(x) is the temperature within region Ω, p(x) is the power density of the heat source, and k(x) is the thermal conductivity of the material. On the first kind of Dirichlet boundary Γ1, the temperature T(x) is a known φ1(x); on the second kind of 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 third type of Robin boundary Γ3, also known as the convective boundary, corresponds to the interface that exchanges heat with the surrounding fluid, such as the surface of a cooling fan and radiator. (Parameters...) Where h(x) is the thermal convection coefficient at the third type Robin boundary, and φ3(x) is a given value.

[0052] In step 2 of this invention, based on the Feynman-Kac formula, a probabilistic solution to the steady-state thermal analysis problem of a three-dimensional chip is obtained using standard reflection Brownian motion.

[0053] Specifically, considering a region that satisfies (1), according to the Feynman-Kac formula, the probabilistic solution to this mixed boundary value problem, i.e., the temperature value T(x) at point x, is [9]:

[0054]

[0055] Among them, E x Let X be the mathematical expectation. t It is a reflection Brownian motion, where L(t) is the corresponding local time. It is the moment when the reflected Brownian motion reaches the Dirichlet boundary of the first kind. It is the position where the reflected Brownian motion reaches the first kind of Dirichlet boundary. The integral terms in (2)-(4) represent the contributions of the right-hand side of the Poisson equation to the temperature of the point to be measured under the first kind of Dirichlet boundary condition, the second kind of Neumann boundary condition, and the third kind of Robin boundary condition, respectively.

[0056] Feynman-Kac functionals Defined as

[16] :

[0057]

[0058] The reflection Brownian motion and its corresponding local time are obtained through the Skorokhod problem based on the definition of Brownian motion (BM). Assume D is a function of R. 3 One of them has C 1 The boundary domain, the Skorokhod problem is defined as follows: Let f∈C([0,∞), R 3 ), is from [0, ∞) to R 3 A continuous function. A pair of (ξ) t L t ) is a solution to the Skorokhod equation S(f;D), if:

[0059] 1.ξ in The middle part is continuous;

[0060] 2. L(t) is a non-decreasing function if and only if Increased over time, that is

[0061]

[0062] 3. The Skorokhod equation holds:

[0063]

[0064] Where n(x) represents Unit outward normal vector at point, index function Defined as:

[0065]

[0066] Suppose f(t) is from The initial Brownian motion path, and (X) t L t If X is a solution to the Skorokhod problem S(f;D), then X t It is a reflective Brownian motion starting from x. The local time is not an independent process, but rather related to the standard reflective Brownian motion X. t Related, defined as

[17] :

[0067]

[0068] Among them, D ∈ It includes The width of the band region is ∈ and L(t) is X t The local time of reflection measures the Brownian motion X. t The vanishing neighborhood D of the second-type Neumann or third-type Robin reflection boundary within the time interval [0,t] ∈ The amount of time spent here

[0069] Therefore, the Feynman-Kac formula provides a method for utilizing the reflection of Brownian motion (RBM) X. t Methods for solving mixed boundary value problems of the Poisson equations, particularly local probabilistic solutions involving Neumann boundaries of the second kind or Robin boundaries of the third kind, where the probabilistic solution is given by the reflection Brownian motion X. tThe expectation of the path integral is calculated. Among them, (2)-(4) are the integral terms on the path with respect to the first kind of Dirichlet boundary conditions, the second kind of Neumann boundary conditions, the third kind of Robin boundary conditions, and the right-hand side terms of the Poisson equation, respectively.

[0070] Typically, reflective Brownian motion is simulated using the WOS method. The starting point x of the reflective Brownian motion is the center of the sphere. According to the Green's function of the sphere, the position of the next point in the reflective Brownian motion is uniformly distributed on the sphere's surface. Therefore, a series of walks or jumps on the sphere can be used to simulate the path of the reflective Brownian motion. To speed up the WOS process, the sphere radius for each jump is set to its maximum value, i.e., the distance from the current point to the boundary. When a point in the reflective Brownian motion collides with a first-type Dirichlet boundary, it is absorbed by the boundary, and the reflective Brownian motion terminates. When a point in the reflective Brownian motion collides with a second-type Neumann boundary or a third-type Robin boundary, it is reflected inward by the boundary. For details on the calculation of the local time, see step 4.3.

[0071] In step 3 of this invention, all input parameters are transmitted to the GPU, N GPU threads are started according to the required number of sampling paths N, and the program is initialized.

[0072] Specifically, a GPU (Graphics Processing Unit) is a processor designed specifically for parallel computing. Taking NVIDIA GPUs as an example, these GPUs employ an architecture supporting the CUDA (Compute Unified Device Architecture) programming model. Their core advantage lies in having thousands of stream processors, enabling them to execute a large number of threaded tasks simultaneously. Unlike CPUs, which focus on complex logic and low latency, GPUs use a SIMT (Single Instruction Multiple Threads) execution model. By breaking down computational tasks into numerous parallelizable subtasks, they achieve data-level parallelism. For example, when processing image pixels or performing matrix operations, each stream processor can independently process different data units, significantly improving throughput.

[0073] GPU parallel computing is achieved through hierarchical thread organization: threads form thread blocks, and multiple thread blocks form a grid. At the hardware level, Stream Multiprocessors (SMs) dynamically schedule warps to execute the same instruction stream, alleviating memory access bottlenecks through multi-level storage structures such as registers and shared memory. This architecture is particularly suitable for handling computational tasks with high parallelism and low data dependencies.

[0074] When starting the GPU randomized solution program, the total number of threads is specified as the required number of paths, N, with each path corresponding to a specific thread on the GPU. In the i-th thread, several variables are used to store the current information of the path corresponding to that thread, including the current jump step number np and the current Feynman-Kac functional. Value, current position X t When the i-th thread initializes its program, it sets its current jump step number X. t Initialized to 0, the current Feynman-Kac functional The value is initialized to 1, and the current position X0 is initialized to x. Additionally, during program initialization, the number of jump steps k that have been completed in multiples of δ on all paths is initialized to 0. Space is reserved in the GPU's global memory to store the cumulative temperature values ​​contributed by the third type Robin boundary and heat source region when the number of jump steps exceeds kδ. Then initialize the data in it to 0.

[0075] In step 4 of this invention, based on the current point X t At the given location, a corresponding random walk process is executed in the GPU; random walks are performed using WOS and WOG in the passive and hot regions, respectively; for three different types of boundary conditions, direct termination, local time calculation, and local time calculation and update are adopted, respectively. Three different calculation strategies are used, and the temperature value of the current path is accumulated based on the discrete form of the probabilistic solution;

[0076] Specifically, if the current point is in the passive region, a random walk (WOS) is performed, i.e., jump to step 4.1; if the current point is on a Dirichlet boundary of the first type... If the current point is within the area, the random walk terminates, i.e., jump to step 4.2; if the current point is on the second type of Neumann boundary... Within the region, calculate the local time and accumulate the temperature value, i.e., jump to step 4.3; if the current point X t In the third type of Robin boundary Within the region, calculate the local time and update. The temperature value is accumulated, i.e., proceed to step 4.4; if the current point is within the heat source area, a WOG (Walk-Operate-Ground) is used to perform a random walk and accumulate the temperature value, i.e., proceed to step 4.5. All jump rules in steps 4.1 to 4.5 are as follows: Figure 3 As shown.

[0077] Because newer GPU architectures (such as the Volta architecture) support independent scheduling of threads within each thread block, all position determination and temperature accumulation are performed in the same GPU function instead of being split into multiple GPU function modules

[11] , thus avoiding the overhead caused by inter-thread communication via global memory on the GPU, such as Figure 4 As shown.

[0078] Step 4.1 of the present invention: If the current point X t In the passive region, WOS is used for random walks.

[0079] Specifically, if the current point X t In the passive region and not at the first type Dirichlet boundary Regions or second-type Neumann boundaries, third-type Robin boundaries Within the area, WOS is used for random walking, with a jump radius equal to the current point X. t The shortest distance to the three boundaries and heat source regions. Uniform sampling is performed on a sphere with this distance as the radius to obtain the next random point X of Brownian motion. t+1 As the current point X for the next step t And increment the jump count np of the current path by 1.

[0080] Step 4.2 of the present invention: If the current point X t First type Dirichlet boundary If the area is within the specified range, the random walk will terminate.

[0081] Specifically, an extremely thin band-like region was constructed near the first type of Dirichlet boundary. Its width is If the current point X t First type Dirichlet boundary If the point is within the area, it will be absorbed onto the boundary by the first type of Dirichlet boundary, and the random walk on this path will terminate.

[0082] This termination time is recorded as Point X t The point closest to the first type of Dirichlet boundary is denoted as The integral term for the temperature contribution of the random walk termination is:

[0083]

[0084] When the third type of Robin boundary is not considered, the Feynman-Kac functional It can be ignored. However, once a third type of Robin boundary is introduced, At this point, it is necessary to record and update the random jump path. The specific calculation details are explained in step 4.4.

[0085] It should be noted that the fitting technique proposed in this invention does not consider paths terminating at the first type of Dirichlet boundary.

[0086] In step 4.3 of this invention, if the current point X t Second type Neumann boundary Within the region, local time is calculated and temperature values ​​are accumulated.

[0087] Specifically, an extremely thin band-like region was constructed near the second-type Neumann boundary. Its width is If the current point X t Second type Neumann boundary Within the region, it is necessary to calculate the local time corresponding to each point.

[0088] If the current point X t Second type Neumann boundary Within the area, WOS is still used for random walks, with the jump radius fixed as a constant Δx or 2Δx. Each jump results in the next random point X obtained through Brownian motion. t+1 As the current point X for the next step t And increment the current path jump count np by 1. Let the current point X be... t Distance R to the second-type Neumann boundary r When 0 < R r When the value is less than Δx, let the jump radius r = 2Δx; otherwise, let the jump radius r = Δx. Using this fixed-step WOS jump method, the local time can be expressed as

[18]

[19] :

[0089]

[0090] Where Δx is a preset value. Fixed jump step size within the strip region, and Entering and exiting the jump point respectively hour The value, The path is from the start to time t. j Staying in the banded area The cumulative number of steps within; it is worth noting that only X t The jump path is located at Only within the banded area Make a contribution;

[0091] exist Within the region, the jump step size of WOS is Δx or 2Δx, therefore the local time corresponding to the jump is (Δx). 2 / 3 or (2Δx) 2 / 3; Therefore, in implementation, the method of the present invention will It is considered an increasing index value; after each step of WOS, The increment of the value relative to the previous step will be 0, 1, or 4, corresponding to the random point not being in the previous step, respectively. The area, in Within the region and with a jump radius of Δx, or In the region with a jump radius of 2Δx, its integral term for the temperature contribution is:

[0092]

[0093] Similarly, when the third type of Robin boundary is not considered, the functional However, once the Robin boundary is introduced, At this point, in equation (13) The calculation details are explained in step 4.4.

[0094] In step 4.4 of this invention, if the current point X t In the third type of Robin boundary Within the region, calculate the local time and update. And accumulate temperature values.

[0095] Specifically, a very thin band-shaped region was constructed near the third type of Robin boundary. If the current point X t In the third type of Robin boundary Within the area, WOS is used for random walking, and each jump yields the next random point X, which is a Brownian motion. t+1 As the current point X for the next step t And increment the current path jump count np by 1, using the same method as in step 4.3.

[0096] In addition, it is also necessary to Perform an update. This update is related to the path in the region. This relates to collisions with the third type of Robin boundary. Random points are located within... When walking randomly within a region, one often leaves without colliding with the boundary. Area. The method of this invention tracks the entry and exit of random points. Each entry time t of the region j，in and export time t j，out And calculate the local time increment dL(t) for each collision that occurs within this time period. This is to calculate the collision point... corresponding It is necessary to record the WOS path for entering and exiting. Coefficient increment between regions and the number of collisions m between the path and the boundary within the same time period. j ,make

[0097]

[0098] Then, It can be approximated as:

[0099]

[0100] Therefore, the integral term of its contribution to temperature is:

[0101]

[0102] Due to parameters The value only appears on the third type of Robin boundary, so The value is only when It changes over time. When At that time, the random walk path reaches the first type of Dirichlet boundary and is absorbed, while This represents the final decay of the path by the third type of Robin boundary during the random process at the end of the random walk.

[0103] Feynman-Kac functionals It is the attenuation of the path caused by the third type of Robin boundary; that is, for a mixed boundary problem with only the first type of Dirichlet boundary and the second type of Neumann boundary, the functional... It equals 1. However, once the third type of Robin boundary is involved, the functional... It will not be 1, and will act as a decay coefficient affecting all integral terms, not just those related to the third type of Robin boundary conditions.

[0104] Therefore, in the process of calculating the path integral, the function It needs to be dynamically updated as a variable each time the path reaches a third-class Robin boundary. The value will suddenly decrease. Furthermore, as the third type of Robin boundary reflection increases, The value gradually changes from 1 to 0, leading to a decrease in the contribution of subsequent integration terms and the eventual convergence of the path integral. It is worth noting that... It is only related to a given parameter c on the third type Robin boundary, so its value only changes when the path reaches the third type Robin boundary.

[0105] In step 4.5 of this invention, if the current point X t Within the heat source area, WOG is used for random walks and to accumulate temperature values.

[0106] Specifically, if the current point X t Within the heat source area, WOG is used for random movement, with each jump leading to the next random point X. t+1As the current point X for the next step t This increments the current path jump count np by 1 and accumulates the temperature value. To implement WOG (Word of the Heat), the heat source region needs to be discretized into a uniform grid, dividing it into several cuboids of the same size, resulting in a discrete grid for the heat source region. Each grid has a corresponding power density and power value, and the grid node is the center point of the grid. Thermal resistance and thermal conduction can be calculated between adjacent grid nodes. For example, if the length and width of two adjacent grids are both h... x and h y The heights are h zi and h zj The thermal conductivity of the two grids are k i and k j Then R(i,j) is the thermal resistance between two grid nodes, defined as:

[0107]

[0108] Where g(i,j) is the thermal conductivity between two grid nodes.

[0109] The probability P(i,j) of jumping from grid node i to the neighboring grid node j is:

[0110]

[0111] Where degree(j) is the number of neighboring grid nodes of grid node i.

[0112] X t The temperature contribution f(X) of the heat source region obtained at node j t In the WOG method, this can be represented as r(j), specifically in the form:

[0113]

[0114] Where p j Let be the power at grid node j. Its integral term for the temperature contribution is:

[0115]

[0116] A virtual grid is constructed at the boundary between the heat source region and the passive region, with a corresponding power of zero. Two different stochastic methods, WOS and WOG, are then switched between. Current point X t When entering this mesh from a passive region, you will be pulled to the center point of the nearest mesh, and then randomly jump to a nearby mesh node or a passive region using the WOG method.

[0117] In step 5 of this invention, it is determined whether the number of jump steps of the current path exceeds an integer multiple of the number of jump steps δ, i.e., np≥kδ. If so, proceed to step 6; otherwise, proceed to step 4.

[0118] Specifically, after each jump in step 4, the jump count np of the i-th path increases by a certain number of steps. If the jump count np exceeds an integer multiple of the jump step count δ (i.e., np ≥ kδ), then the current cumulative temperature value of this path should be recorded, and the process should proceed to step 6; otherwise, the path continues to jump, and the process returns to step 4.

[0119] In step 6 of this invention, the cumulative temperature T is recorded when the number of jump steps is kδ for the i-th path, calculated using the discrete form of the Feynman-Kac formula. i，kδ (x), if all N paths have completed at least kδ jumps, proceed to step 7; otherwise, proceed to step 4.

[0120] Specifically, according to step 4, and following the discrete form of the Feynman-Kac formula

[13] , the real-time cumulative temperature value of the i-th path starting from point x at the jump step number kδ is:

[0121]

[0122] Where j′ represents each jump point in the path, Indicates the second type of Neumann boundary and the third type of Robin boundary. Jump points in the area, This represents the jump point within the heat source region. Note that when saving this temperature value, the actual number of jump steps np in the path may be slightly greater than kδ, but the calculation only calculates the cumulative temperature value from step 0 to step kδ. This temperature value T... i，kδ (x) will be saved in the GPU global memory. After saving, it is determined whether all paths have completed kδ jumps. If so, the cumulative temperature sampling of all kδ steps has been completed, the average value can be calculated, and the process proceeds to step 7. Otherwise, the average value is not calculated for the time being, the current path continues to jump, and the process proceeds to step 4.

[0123] In step 7 of this invention, the average value of the total contribution of the third type Robin boundary and the heat source region to the cumulative temperature value is calculated when the number of jump steps for N paths is kδ. After adding this average value, the temperature is fitted within the temperature fitting window w. To determine whether the fitted temperature has converged, i.e., to determine... If yes, terminate the program and calculate the temperature T(x) at the point x to be determined as the output; otherwise, increment k by 1 and proceed to step 4.

[0124] Specifically, in practical 3D chip thermal analysis, there are typically numerous Neumann boundaries of type II and Robin boundaries of type III. For these two types of boundaries, due to the Feynman-Kac functional... Each path requires a large number of jump steps to converge the accumulated temperature, resulting in a long program runtime. However, in practical thermal analysis applications, when the number of jump steps is large enough, the temperature obtained by PIRW usually begins to converge exponentially. This characteristic can be used to reduce the number of steps required to achieve the desired accuracy through fitting.

[0125] Since the path typically terminates quickly when a first-type Dirichlet absorbing boundary exists in the problem, without affecting the execution time of the parallel program, this invention does not consider the first-type absorbing boundary in the fitting. Furthermore, in practical chip thermal analysis problems, the second-type Neumann boundary condition is usually an adiabatic boundary condition, contributing zero to the cumulative temperature value, and is not considered in the fitting. Therefore, the average of the total contribution of the third-type Robin boundary and the heat source region to the cumulative temperature value for N paths at kδ steps can be written as:

[0126]

[0127] Where I(.) is the indicator function.

[0128] Approximating formula (22), we can obtain:

[0129]

[0130]

[0131] According to the properties of Brownian motion, within a finite space, after a sufficient number of jumps, the spatial distribution of the current positions of all paths will tend to stabilize. Because... It only depends on the spatial distribution of the current positions of all paths at the current step j. After a sufficient number of jump steps, It will stabilize to a constant. Furthermore, due to the decay of the Feynman-Kac functional caused by the exponent... and dL(t j ) control, and and dL(t j The temperature increment is only related to the spatial distribution of the current positions of all paths at the current step j. Therefore, the Feynman-Kac functional decays exponentially with the number of jump steps. Combining these two points, when the number of jump steps is sufficiently large, the temperature increment at each step should be approximately proportional to the current Feynman-Kac functional, i.e., decaying exponentially with the number of jump steps. This leads to the cumulative temperature value beginning to converge exponentially, as shown below. Figure 5 As shown.

[0132] Based on the above assumptions, when all N paths have completed kδ steps of transitions, the average of the total contributions of the third-type Robin boundary and heat source region to the cumulative temperature value at multiple transition step points within the temperature fitting window w is used. The convergence temperature T is estimated by fitting the equation (25). c And parameters a, b, and c. The fitted T c Recorded as

[0133]

[0134] In actual operation, the decision to terminate the program is made by determining whether the fitted values ​​have converged. The current fitted temperature is then set. Compared with the previous fitted temperature If a comparison is made, Then the latest fitted value is used as the final result and all paths are terminated.

[0135] If the fitted values ​​have stabilized, i.e. The program will then terminate and output the temperature at point x.

[0136]

[0137] in, The termination time, Termination point X t The point closest to the first type of Dirichlet boundary, j′ represents each jump point in the path. This indicates a jump point within the second type of reflection boundary region.

[0138] Otherwise, that is Increment k by 1, then proceed to step 4 to continue the random walk.

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

[0140] (1) A parallel high-precision stochastic method for steady-state thermal analysis of three-dimensional chips was implemented on a GPU. The parallel computation of thousands of computing cores on the GPU significantly reduced the computation time of the PIRW program.

[0141] (2) A fitting technique is proposed, which utilizes the property that the temperature value converges approximately exponentially to reduce the number of steps required for the path to converge to temperature in the random method.

[0142] This application addresses steady-state thermal analysis problems under mixed boundary conditions in 3D chips. The proposed GPU parallel acceleration method is a powerful tool for solving mixed boundary value problems of the Laplace and Poisson equations, offering a significant speed advantage for solving complex thermal problems in large-scale 3D chips. Furthermore, this method can be extended to field analysis problems of the Laplace and Poisson equations with mixed boundary conditions in other engineering fields.

[0143] The advantages of this method are:

[0144] (1) Implement a parallel, high-precision stochastic three-dimensional chip steady-state thermal analysis method on GPU, utilizing thousands of computing cores on GPU for parallel computation, significantly reducing the computation time of PIRW program.

[0145] (2) Taking advantage of the property that the temperature value converges approximately exponentially in random walk, a fitting technique is proposed to achieve early termination of random walk, thereby reducing the number of steps required for each path in the random method to reach temperature convergence, so as to achieve the purpose of acceleration. Attached Figure Description

[0146] Figure 1 The algorithm flowchart proposed in this application;

[0147] Figure 2 A schematic diagram of a 3D heat conduction problem with mixed boundary conditions;

[0148] Figure 3 A schematic diagram of path jumps in the PIRW algorithm;

[0149] Figure 4 A diagram illustrating the execution of branch instructions in CPU and GPU environments;

[0150] Figure 5 A schematic diagram illustrating the principle of temperature fitting;

[0151] Figure 6 Test case diagram;

[0152] Figure 7 The results of fitting the temperature. Detailed Implementation

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

[0154] Experimental conditions settings:

[0155] like Figure 6The test case structure consists of a three-layer thermal conductivity stack. The gray areas represent the heat source region that continuously generates heat, which is the transistor layer or memory layer in the chip. The top and bottom cuboid regions are passive regions, simulating the heat sink, dielectric layer, PCB, etc., in the packaged chip. Testing was conducted under three examples. In all examples, the top and bottom surfaces were Robin boundary conditions, and the ambient temperature was T. amb =20℃; four sides are adiabatic boundaries; the heat source area uses a generated 4-core chip example with a total power of 356W, such as Figure 6 As shown; the convection coefficient at the bottom Robin boundary is fixed at h = 8700 W / m. 2 ·℃. The remaining parameters are shown in the table below:

[0156]

[0157] Table 1: Parameter Settings for the Three-Layer Stacked Region

[0158] In this method, the number of paths MN = 2.5e4, Δx = 5e-7 meters. ε th =1e-4, jump step interval δ = 4e5 steps, temperature fitting window w = 4. For the PIRW simulator used for comparison without GPU and fitting acceleration, the convergence condition is jump step number NP ≥ 1e7 steps or Feynman-Kac functional. The reference temperature was obtained by high-precision discretization using the commercial finite element analysis software ANSYS, and the solution speed was referenced to PIRW, which was executed on the CPU and did not use fitting techniques.

[0159] The parameters used for fitting, including the minimum, maximum, and initial values, are shown in the table below:

[0160] Table 2: Parameter settings used for fitting

[0161]

[0162] Implementation Example 1

[0163] The objective of this example is to compare the accuracy of the accelerated solution method proposed in this application with that of commercial software, as well as the improvement in solution speed compared to the original unaccelerated method.

[0164] For each example above, two points were selected for simulation: one was a hotspot, and the other was a randomly selected point. Each point was run 10 times, and the average and standard deviation of the cumulative temperature values, as well as the average number of steps and the time taken, were recorded.

[0165] The accuracy and speed results of hotspot temperatures obtained using the acceleration method of this application are shown in Table 3.

[0166] Regarding accuracy comparisons, across all six test points, the temperature error of the unaccelerated CPU-based PIRW emulator ranged from 0.23K to 1.86K, while the temperature error of the GPU-accelerated PIRW emulator ranged from 0.16K to 1.75K. Therefore, GPU acceleration did not result in any loss of accuracy for the PIRW emulator.

[0167] Furthermore, the temperature error of the PIRW emulator using fitting techniques and GPU parallel acceleration is 0.17K-0.89K, which is smaller than the temperature errors of the CPU-based PIRW emulator and the GPU-based PIRW emulator without fitting. Therefore, the fitting technique does not introduce additional error.

[0168] In some cases, fitting techniques can even reduce errors. This is because when the geometry becomes larger or the convection coefficient of a third-type Robin boundary decreases, the decay of the Feynman-Kac functional slows down, increasing the number of jump steps required for PIRW to converge. This could exceed the maximum jump step limit, causing the program to stop prematurely before convergence, resulting in a larger error in the accumulated temperature value. Using fitting techniques, however, the program may stop and obtain the fitted temperature value with fewer jumps.

[0169] Regarding speed comparison, for all six test points, the acceleration method proposed in this application, combining GPU acceleration and fitting techniques, achieves an average speedup of 64.1 times compared to the original PIRW program. Specifically, using GPU acceleration alone achieves an average speedup of 25.3 times; combining fitting techniques, while maintaining similar accuracy to GPU acceleration without fitting, further increases the average speedup by 2.53 times by reducing the required number of jump steps. Therefore, the method in this application can effectively accelerate high-precision stochastic three-dimensional chip steady-state thermal analysis.

[0170] Table 3: Comparison of speed and accuracy between the accelerated method in this application and the original method

[0171]

[0172] Implementation Example 2

[0173] The goal of this example is to verify the effectiveness of the fitting technique.

[0174] In the execution of case 1-1, this application uniformly sets the number of path steps to 1e7 and records the cumulative temperature values ​​obtained with and without using fitting techniques at certain jump points, such as... Figure 7As shown in the figure, it can be observed that the temperature values ​​converge much faster after using the fitting technique. The fitting result obtained in approximately 2e6 steps is already close to the convergence value obtained in approximately 1e7 steps without using the fitting result. Therefore, the experiment verifies that the fitting technique can reduce the number of steps required for the cumulative temperature values ​​to converge, thereby further accelerating the operation of the PIRW simulator.

Claims

1. A high-precision stochastic method for accelerating steady-state thermal analysis of three-dimensional chips based on GPU parallelism and path fitting techniques, characterized in that, It includes, First, the steady-state thermal analysis problem of the three-dimensional chip is transformed into a boundary value problem of the Poisson equation with mixed boundary conditions; then, based on the Feynman-Kac formula, the probabilistic solution of the steady-state thermal analysis problem of the three-dimensional chip is obtained by using standard reflection Brownian motion. Third, transmit all input parameters to the GPU, start N GPU threads according to the required number of sampling paths N, and initialize the program; Fourth, based on the current point X t At the given location, a corresponding random walk process is executed in the GPU; random walks are performed using WOS and WOG in the passive and hot regions, respectively; for three different types of boundary conditions, direct termination, local time calculation, and local time calculation and update are adopted, respectively. Three different calculation strategies; Fifth, determine whether the number of jump steps in the current path exceeds an integer multiple of the jump step number δ, i.e., np≥kδ. If so, proceed to step six; otherwise, proceed to step four; Sixth, record the cumulative temperature T of the i-th path when the number of jump steps is kδ, calculated using the discrete form of the Feynman-Kac formula. i,kδ (x), if all N paths have completed at least kδ steps of jumps, proceed to step 7; otherwise, proceed to step 4. Seventh, calculate the average total contribution of the third-type Robin boundary and heat source region to the cumulative temperature value when the number of jump steps for N paths is kδ. After adding this average value, the temperature is fitted within the temperature fitting window w. To determine whether the fitted temperature has converged, i.e., to determine... If yes, terminate the program and calculate the temperature T(x) at the point x to be determined as the output; otherwise, increment k by 1 and proceed to step four. The specific methods include: Input parameters: region Ω and its boundaries Geometric information, region boundaries The three types of boundaries Γ i and boundary conditions i = 1, 2, 3, near the first-kind absorbing boundary Thickness of the strip region Near the boundaries of Type II and Type III reflections Thickness of the strip region and Within the strip region, the WOS jump radius Δx, the location of the point x to be determined, the number of sampling paths N, the interval jump steps δ for initiating fitting, and the convergence criterion ε for judging the fitted value. th Temperature fitting window w; Output: Temperature T(x) at point x; Including specific steps: Step 1: Transform the steady-state thermal analysis problem of the 3D chip into a boundary value problem of the Poisson equation with mixed boundary conditions; Step 2: Based on the Feynman-Kac formula, and using standard reflection Brownian motion, obtain the probabilistic solution to the three-dimensional chip steady-state thermal analysis problem; Step 3: Transfer all input parameters to the GPU, start N GPU threads according to the required number of sampling paths N, and initialize the program; Step 4: Based on the current point X t At the given location, a corresponding random walk process is executed in the GPU; random walks are performed using WOS and WOG in the passive and hot regions, respectively; for three different types of boundary conditions, direct termination, local time calculation, and local time calculation and update are adopted, respectively. Three different computational strategies; Step 5: Determine whether the number of jump steps in the current path exceeds an integer multiple of the number of jump steps δ, i.e., np≥kδ. If yes, proceed to step 6; otherwise, proceed to step 4. Step 6: Record the cumulative temperature T for the i-th path when the number of jump steps is kδ, calculated using the discrete form of the Feynman-Kac formula. i，kδ (x), if all N paths have completed at least kδ jumps, proceed to step 7; otherwise, proceed to step 4. Step 7: Calculate the average total contribution of the third-type Robin boundary and heat source region to the cumulative temperature value when the number of jump steps for N paths is kδ. After adding this average value, the temperature is fitted within the temperature fitting window w. To determine whether the fitted temperature has converged, i.e., to determine... If yes, terminate the program and calculate the temperature T(x) at the point x to be determined as the output; otherwise, increment k by 1 and proceed to step 4.

2. The method according to claim 1, characterized in that, In step 1, the steady-state thermal analysis problem of the three-dimensional chip is transformed into a boundary value problem of the Poisson equation with mixed boundary conditions; The problem to be solved is steady-state thermal analysis of a region containing an internal heat source; the three-dimensional region for steady-state thermal analysis is... Within region Ω, the temperature field T(x) satisfies the Poisson equation, and the region boundary... The boundary is a mixture, meaning it contains a first-kind Dirichlet boundary Γ1, a second-kind Neumann boundary Γ2, and a third-kind Robin boundary Γ3. Given these three types of boundary conditions, the field differential equations and boundary conditions are as follows: Where T(x) is the temperature within the region Ω, p(x) is the power density of the heat source, 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 third type of Robin boundary Γ3, also known as the convective boundary, corresponds to the interface that exchanges heat with the surrounding fluid, such as the surface of a cooling fan and radiator. (Parameters...) Where h(x) is the thermal convection coefficient at the third type Robin boundary, and φ3(x) is a given value.

3. The method according to claim 1, characterized in that, In step 2, based on the Feynman-Kac formula, the probabilistic solution to the three-dimensional chip steady-state thermal analysis problem is obtained using standard reflection Brownian motion. Considering a region satisfying (1), according to the Feynman-Kac formula, the probabilistic solution to this mixed boundary value problem, i.e., the temperature value T(x) at point x, is: Among them, E x Let X be the mathematical expectation. t It is a reflection Brownian motion, where L(t) is the corresponding local time. It is the moment when the reflected Brownian motion reaches the Dirichlet boundary of the first kind. It is the position where the reflected Brownian motion reaches the first kind of Dirichlet boundary; the integral terms in (2)-(4) respectively represent the contributions of the right-hand side of the Poisson equation to the temperature of the point to be measured under the first kind of Dirichlet boundary condition, the second kind of Neumann boundary condition, and the third kind of Robin boundary condition; Feynman-Kac functionals Defined as: The Brownian motion of reflection and the corresponding local time are obtained from the Skorokhod problem based on the definition of Brownian motion; assuming D is R 3 One of them has C 1 The boundary domain, the Skorokhod problem is defined as follows: Let f∈C([0,∞), R 3 ), is from [0, ∞) to R 3 A continuous function; a pair of (ξ) t L t ) is a solution to the Skorokhod equation S(f;D), if: 1). ξ in The middle part is continuous; 2) L(t) is a non-decreasing function if and only if Increased over time, that is 3). The Skorokhod equation holds: Where n(x) represents Unit outward normal vector at point, index function Defined as: Suppose f(t) is from The initial Brownian motion path, and (X) t L t If X is a solution to the Skorokhod problem S(f;D), then X t It is a reflective Brownian motion starting from x; the local time is not an independent process, but rather related to the standard reflective Brownian motion X. t Related, defined as: Among them, D ∈ It includes The width of the band region is ∈ and L(t) is X t The local time of reflection measures the Brownian motion X. t The vanishing neighborhood D of the second-type Neumann or third-type Robin reflection boundary within the time interval [0, t] ∈ The amount of time spent here The Feynman-Kac formula provides a method for utilizing the Brownian motion of reflection X. t Methods for solving mixed boundary value problems of the Poisson equations, particularly local probabilistic solutions involving Neumann boundaries of the second kind or Robin boundaries of the third kind, where the probabilistic solution is given by the reflection Brownian motion X. t The expectation of the path integral is calculated; where (2)-(4) are the integral terms on the path with respect to the first type of Dirichlet boundary conditions, the second type of Neumann boundary conditions, the third type of Robin boundary conditions, and the right-hand side terms of the Poisson equation, respectively. The reflection Brownian motion is simulated using the WOS method. The starting point x of the reflection Brownian motion is the center of the sphere. According to the Green's function of the sphere, the position of the next point in the reflection Brownian motion is uniformly distributed on the surface of the sphere. Therefore, a series of walking or jumping on the sphere can be used to simulate the path of the reflection Brownian motion. At the same time, in order to speed up the WOS process, the radius of the sphere for each jump is set to the maximum value, that is, the distance from the current point to the boundary. When the point of the reflection Brownian motion collides with the first type of Dirichlet boundary, the point will be absorbed by the boundary, and the reflection Brownian motion will terminate. When the point of the reflection Brownian motion collides with the second type of Neumann boundary or the third type of Robin boundary, the point will be reflected inward by the boundary. The specific details of the local time calculation are as described in step 4.

3.

4. The method according to claim 1, characterized in that, In step 3, all input parameters are transmitted to the GPU, N GPU threads are started according to the required number of sampling paths N, and the program is initialized. Specifically, GPUs are processors designed specifically for parallel computing. Taking NVIDIA GPUs as an example, these GPUs use an architecture that supports the CUDA programming model. Their core advantage lies in having thousands of stream processors, enabling them to execute a large number of threaded tasks simultaneously. Unlike CPUs, which focus on complex logic and low latency, GPUs use a Single Instruction Multithreaded (SIMT) execution model. By breaking down computational tasks into a large number of parallelizable subtasks, they achieve data-level parallelism. For example, when processing image pixels or performing matrix operations, each stream processor can independently process different data units, significantly improving throughput. GPU parallel computing is achieved through hierarchical thread organization: threads form thread blocks, and multiple thread blocks form a grid; at the hardware level, the Streaming Multiprocessor (SM) dynamically schedules thread bundles (warps) to execute the same instruction stream, and uses multi-level storage structures such as registers and shared memory to alleviate memory access bottlenecks; this architecture is suitable for handling computing tasks with high parallelism and low data dependency. When starting the GPU randomized solution program, the total number of threads is specified as the required number of paths N, and each path corresponds to a specific thread on the GPU. In the i-th thread, several variables are used to store the current information of the path corresponding to that thread, including the current jump step np and the current Feynman-Kac functional. Value, current position X t When the i-th thread initializes its program, it sets its current jump step number X. t Initialized to 0, the current Feynman-Kac functional The value is initialized to 1, and the current position X0 is initialized to x; in addition, during program initialization, the number of jump steps k that have been completed in multiples of δ for all paths is initialized to 0, and space is reserved in the GPU global memory to store the cumulative temperature value contributed by the third type Robin boundary and heat source region when the number of jump steps exceeds kδ. Then initialize the data in it to 0.

5. The method according to claim 1, characterized in that, In step 4, based on the current point X t At the given location, a corresponding random walk process is executed in the GPU; random walks are performed using WOS and WOG in the passive and hot regions, respectively; for three different types of boundary conditions, direct termination, local time calculation, and local time calculation and update are adopted, respectively. Three different calculation strategies are used, and the temperature value of the current path is accumulated according to the discrete form of the probabilistic solution; Specifically, if the current point is in a passive region, a random walk (WOS) is performed, i.e., jump to step 4.1; if the current point is on a first-type Dirichlet boundary... If the current point is within the area, the random walk terminates, i.e., jump to step 4.2; if the current point is on the second type of Neumann boundary... Within the region, calculate the local time and accumulate the temperature value, i.e., jump to step 4.3; if the current point X t In the third type of Robin boundary Within the region, calculate the local time and update. Accumulate the temperature value, i.e., jump to step 4.4; if the current point is within the heat source area, use WOG to perform random walk and accumulate the temperature value, i.e., jump to step 4.5; Because newer GPU architectures support independent scheduling of threads within each thread block, all position determination and temperature accumulation are performed in the same GPU function instead of being split into multiple GPU function modules, thus avoiding the overhead caused by inter-thread communication through global memory on the GPU. In step 4.1, if the current point X t In the passive region, WOS is used for random walks; Including, if the current point X t In the passive region and not at the first type Dirichlet boundary Regions or second-type Neumann boundaries, third-type Robin boundaries Within the area, WOS is used for random walking, with a jump radius equal to the current point X. t Find the shortest distance to the three boundaries and the heat source region; uniformly sample on a sphere with this distance as the radius to obtain the next random point X of Brownian motion. t+1 As the current point X for the next step t And increment the jump count np of the current path by 1; In step 4.2, if the current point X t First type Dirichlet boundary Within the area, the random walk terminates; Specifically, an extremely thin band-like region was constructed near the first type of Dirichlet boundary. Its width is If the current point X t First type Dirichlet boundary If the point is within the area, it will be absorbed onto the boundary by the first type of Dirichlet boundary, and the random walk on this path will terminate. This termination time is recorded as Point X t The point closest to the first type of Dirichlet boundary is denoted as The integral term for the temperature contribution of the random walk termination is: When the third type of Robin boundary is not considered, the Feynman-Kac functional It can be ignored; however, once a third type of Robin boundary is introduced, At this point, it is necessary to record and update the random jump path. The specific calculation details are explained in step 4.4; It should be noted that the fitting technique proposed in this invention does not consider paths terminating at the first type of Dirichlet boundary; In step 4.3, if the current point X t Second type Neumann boundary Within the region, local time is calculated and temperature values ​​are accumulated; Specifically, an extremely thin band-like region was constructed near the second-type Neumann boundary. Its width is If the current point X t Second type Neumann boundary Within the region, it is necessary to calculate the local time corresponding to each point; If the current point X t Second type Neumann boundary Within the area, WOS is still used for random walks, with the jump radius fixed as a constant Δx or 2Δx. Each jump results in the next random point X obtained through Brownian motion. t+1 As the current point X for the next step t And increment the current path jump count np by 1; let the current point X be... t Distance R to the second-type Neumann boundary r When 0 < Rr < Δx, let the jump radius r = 2Δx; otherwise, let the jump radius r = Δx. Using this fixed-step WOS jump method, the local time can be expressed as: Where Δx is a preset value. Fixed jump step size within the strip region, and Entering and exiting the jump point respectively hour The value, The path is from the start to time t. j Staying in the banded area The cumulative number of steps within; it is worth noting that only X t The jump path is located at Only within the banded area Make a contribution; exist Within the region, the jump step size of WOS is Δx or 2Δx, therefore the local time corresponding to the jump is (Δx). 2 / 3 or (2Δx) 2 / 3; Therefore, in implementation, this method will It is considered an increasing index value; after each step of WOS, The increment of the value relative to the previous step will be 0, 1, or 4, corresponding to the random point not being in the previous step, respectively. The area, in In the region and the jump radius is Δx, or In the region with a jump radius of 2Δx, its integral term for the temperature contribution is: Similarly, when the third type of Robin boundary is not considered, the functional However, once the Robin boundary is introduced, At this point, in equation (13) The calculation details are explained in step 4.4; In step 4.4, if the current point X t In the third type of Robin boundary Within the region, calculate the local time and update. And accumulate temperature values; Specifically, a very thin band-shaped region was constructed near the third type of Robin boundary. If the current point X t In the third type of Robin boundary Within the area, WOS is used for random walking, and each jump yields the next random point X, which is a Brownian motion. t+1 As the current point X for the next step t And increment the current path jump count np by 1, using the same method as in step 4.3; In addition, it is also necessary to Perform an update; this update is related to the path in the region. The collision situation is related to the boundary of the third type of Robin; random points are in When walking randomly within a region, one often leaves without colliding with the boundary. Area; This method tracks the entry and exit of random points. Each entry time t of the region j,in and export time t j，out And calculate the local time increment dL(t) corresponding to each collision that occurs within this time period; in order to calculate the collision point corresponding It is necessary to record the WOS path for entering and exiting. Coefficient increment between regions and the number of collisions m between the path and the boundary within the same time period. j ,make Then, It can be approximated as: Therefore, the integral term of its contribution to temperature is: Due to parameters The value only appears on the third type of Robin boundary, so The value is only when Changes occur over time; when At that time, the random walk path reaches the first-type Dirichlet boundary and is absorbed, while This represents the final decay of the third type Robin boundary on the path during the random process at the end of the random walk; Feynman-Kac functionals It is the attenuation of the path caused by the third type of Robin boundary, that is, for the mixed boundary problem with only the first type of Dirichlet boundary and the second type of Neumann boundary, the functional It equals 1; however, once the third type of Robin boundary is involved, the functional... It will not be 1, and will act as a decay coefficient affecting all integral terms, not just those related to the third type of Robin boundary conditions; Therefore, in the process of calculating the path integral, the function It needs to be dynamically updated as a variable; each time the path reaches the third type Robin boundary, The value will suddenly decrease; and, with the increase of the third type of Robin boundary reflection, The value gradually changes from 1 to 0, leading to a decrease in the contribution of subsequent integration terms and the eventual convergence of the path integral; it is worth noting that, It is only related to the given parameter c on the third type Robin boundary, so its value only changes when the path reaches the third type Robin boundary; In step 4.5, if the current point X t Within the heat source area, WOG is used for random walks and temperature values ​​are accumulated; Specifically, if the current point X t Within the heat source area, WOG is used for random movement, with each jump leading to the next random point X. t+1 As the current point X for the next step t This increments the current path jump count np by 1 and accumulates the temperature value. To implement WOG (Word of the Heat), the heat source region needs to be discretized into a uniform grid, dividing it into several cuboids of the same size, resulting in a discrete grid for the heat source region. Each grid has a corresponding power density and power value, and the grid node is the center point of the grid. Thermal resistance and thermal conduction can be calculated between two adjacent grid nodes. For example, if the length and width of two adjacent grids are both h... x and h y The heights are h zi and h zj The thermal conductivity of the two grids are k i and k j Then R(i,j) is the thermal resistance between two grid nodes, defined as: Where g(i,j) is the thermal conductivity between two grid nodes; The probability P(i,j) of jumping from grid node i to the neighboring grid node j is: Where degree(j) is the number of neighboring grid nodes of grid node i; X t The temperature contribution f(X) of the heat source region obtained at node j t In the WOG method, this can be represented as r(j), specifically in the form: Where p j Let be the power at grid node j; its integral term for the temperature contribution is: A virtual grid is constructed at the boundary between the heat source region and the passive region, with a corresponding power of zero, and two different stochastic methods, WOS and WOG, are switched; the current point X t When entering this mesh from a passive region, you will be pulled to the center point of the nearest mesh, and then randomly jump to a nearby mesh node or a passive region using the WOG method.

6. The method according to claim 1, characterized in that, In step 6, the cumulative temperature T is recorded when the number of jump steps is kδ for the i-th path, calculated using the discrete form of the Feynman-Kac formula. i，kδ (x); Specifically, according to step 4, and following the discrete form of the Feynman-Kac formula, the real-time cumulative temperature value of the i-th path starting from point x at the time of kδ jump steps is: Where j' represents each jump point in the path, Indicates the second type of Neumann boundary and the third type of Robin boundary. Jump points in the area, This represents the jump point within the heat source region; note that when saving this temperature value, the actual number of jump steps np in the path may be slightly greater than kδ, but the calculation only calculates the cumulative temperature value from step 0 to step kδ; this temperature value T i，kδ (x) will be saved in the GPU global memory; after saving, it is determined whether all paths have completed kδ jumps. If so, the cumulative temperature sampling of all kδ steps has been completed, the average value can be calculated, and the process proceeds to step 7. Otherwise, the average value is not calculated for the time being, the current path continues to jump, and the process proceeds to step 4.

7. The method according to claim 1, characterized in that, In step 7, the average value of the total contribution of the third-type Robin boundary and the heat source region to the cumulative temperature value is calculated when the number of jump steps for N paths is kδ. After adding this average value, the temperature is fitted within the temperature fitting window w. To determine whether the fitted temperature has converged, i.e., to determine... If yes, terminate the program and calculate the temperature T(x) at the point x to be determined as the output; otherwise, increment k by 1 and go to step 4. Specifically, in actual 3D chip thermal analysis, there are usually a large number of second-type Neumann boundaries and third-type Robin boundaries; for these two types of boundaries, due to the Feynman-Kac functional... Each path requires a large number of jump steps to bring the accumulated temperature to converge, which results in a long program runtime. However, in practical thermal analysis applications, when the number of jump steps is large enough, the temperature obtained by PIRW usually begins to converge exponentially. This characteristic can be reduced by fitting to reduce the number of steps required to achieve the desired accuracy. Since the path typically terminates quickly when a first-type Dirichlet absorbing boundary exists in the problem, without affecting the execution time of the parallel program, this method does not consider the first-type absorbing boundary in the fitting. Furthermore, in practical chip thermal analysis problems, the second-type Neumann boundary condition is usually an adiabatic boundary condition, contributing zero to the cumulative temperature value, and is not considered in the fitting. Therefore, the average of the total contribution of the third-type Robin boundary and the heat source region to the cumulative temperature value for N paths at kδ steps can be written as: Where I(.) is the indicator function; Approximating formula (22), we can obtain: According to the properties of Brownian motion, within a finite space, after a sufficient number of jumps, the spatial distribution of the current positions of all paths will tend to stabilize; because It only depends on the spatial distribution of the current positions of all paths at the current step j. After a sufficient number of jump steps, It will stabilize to a constant; and because the decay of the Feynman-Kac functional is due to the exponential... and dL(t j ) control, and and dL(t j The temperature increment is only related to the spatial distribution of the current position of all paths at the current step j. Therefore, the Feynman-Kac functional decays exponentially with the increase of the number of jump steps. Based on the above two points, when the number of jump steps is large enough, the temperature increment of each step should be approximately proportional to the current Feynman-Kac functional, that is, it decays exponentially with the increase of the number of jump steps, which leads to the cumulative temperature value starting to converge exponentially. Based on the above assumptions, when all N paths have completed kδ steps of transitions, the average of the total contributions of the third-type Robin boundary and heat source region to the cumulative temperature value at multiple transition step points within the temperature fitting window w is used. i = 0, 1, ..., w-1, the convergence temperature T is estimated by fitting the equation (25). c And parameters a, b, and c; the fitted T c Recorded as In actual operation, the decision to terminate the program is made by determining whether the fitted values ​​have converged; the current fitted temperature is then set. Compared with the previous fitted temperature If a comparison is made, Then the latest fitted value is used as the final result and all paths are terminated; If the fitted values ​​have stabilized, i.e. The program will then terminate and output the temperature at point x. in, The termination time, Termination point X t The point closest to the first type of Dirichlet boundary, j′ represents each jump point in the path. Indicates the jump point within the second type of reflection boundary region; Otherwise, that is Increment k by 1, then proceed to step 4 to continue the random walk.

Images