一种基于GPU并行和拟合的高精度随机法三维芯片稳态热分析加速方法

本发明属于集成电路计算机辅助设计/电子设计自动化(computer aideddesign/electronic designautomatic cad/eda)，具体涉及三维芯片中混合边界条件下稳态热分析问题求解。本方法利用gpu并行和拟合技术，能有效加速处理有关robin边界相关的feynman-kac泛函衰减较慢的难题，进一步提升有关并行效率。

背景技术：

1、现有技术中，随着特征尺寸不断缩小和集成电路设计发展，chiplet/3d ic成为重要技术方向之一，但其中的热问题对集成电路长期可靠性和性能有显著影响。因此，准确高效的热分析变得不可或缺，以确保芯片设计和芯片运行符合温度约束。与此同时，gpu由于其大量并行计算资源而被广泛使用，已被证实在大量电子设计自动化(eda)问题中非常高效。

2、现有的热仿真方法可分为三类：确定性方法、格林函数法、基于神经网络的方法和随机方法。

3、确定性方法包括有限元方法(fem)[1]、[2]和有限差分方法(fdm)[3]、[4]。有限元方法通过将连续区域离散为相互连接的单元，利用变分原理求解偏微分方程。有限差分方法则基于网格点处的微分近似，通过差商代替导数构建代数方程组。虽然确定性方法精度高且适应性强，但求解大规模问题，由于一般涉及大规模线性方程组求解而非常耗时。

4、格林函数法[5]、[6]将场的计算转化为格林函数和源的卷积，但格林函数在实际问题中通常没有解析形式，数值解需要通过确定性方法获得。

5、基于神经网络的方法直接拟合精确标签[7]或使用物理信息损失函数[8]，通过网络预测得到任意点的估计温度值。然而，它们难以处理几何变化的问题，而且网络训练时间长、无法估计误差范围。

6、随机方法作为一种高精度热点分析的潜在解决方案[12]、[13]、[15]，已逐渐在芯片热仿真领域展现其独特优势。基于feynman-kac公式[9]，随机法通过从待求点释放许多粒子进行布朗运动，并计算其路径积分的期望值，获得偏微分方程的单点解。随机法已广泛用于电容提取[10]，并在商业软件中得到广泛实施。路径积分随机行走(pirw)[13]是一种高精度的基于随机方法的求解器，它使用严格准确的局部时间计算和feynman-kac泛函[14]、[16]处理混合边界条件。

7、然而，由于随机法求解结果一般具有的方差，其中n是路径数量[12]，因此需要超过数千条独立路径才能达到高精度。若能利用并行技术，则可大幅提升程序计算效率。然而，由于第三类robin边界在芯片热分析问题中普遍存在，通常每条随机行走路径需要执行数百万次跳跃才能确保所有路径最终温度收敛。这一特性对pirw程序的并行效果带来限制，即便所有路径全部并行，最终并行执行时间还是无法有效缩减。

8、为了提高pirw程序的并行性能，本发明拟提供一种在温度尚未完全收敛时提前终止随机行走，并通过拟合估计最终收敛温度的方法，本方法能进一步提升并行效率，确定提前终止随机行走路径的时机。

9、本方法针对三维芯片中混合边界条件下稳态热分析问题，提出的gpu并行加速方法是求解laplace方程和poisson方程混合边值问题的有力工具，对于解决大规模三维芯片复杂热问题具有显著的速度优势。并且，本方法可以拓展到其他工程领域中具有混合边界条件的laplace方程和poisson方程的场分析问题中。

10、于本发明相关的参考文献；

11、[1]incansys et al.ansys.https://www.ansys.com/,2018.

12、[2]inc.comsol.comsol.https://www.comsol.com/,2022.

13、[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.

14、[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.

15、[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.

16、[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.

17、[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.

18、[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.

19、[9]ean-paulmorillon.numerical solutions of linear mixed boundaryvalue problems using stochastic representations.international journal fornumerical methods in engineering,40(3):387–405,1997.

20、[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.

21、[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.

22、[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.

23、[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.

24、[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.

25、[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.

26、[16]vassilis g papanicolaou.1990.the probabilistic solution of thethird boundary value problem for second order elliptic equations.probabilitytheory and related fields 87,1(1990),27–77.

27、[17]pei hsu.1984.reflecting brownian motion,boundary local time andthe neumann problem.ph.d.dissertation.stanford university.

28、[18]yijing zhou and wei cai.2016.numerical solution of the robinproblem of laplace equations with a feynman–kac formula and reflectingbrownian motions.journal ofscientific computing 69,1(2016),107–121.

29、[19]yijing zhou,wei cai,and elton hsu.2017.computation of the localtime of reflecting brownian motion and the probabilistic representation ofthe neumann problem.communications in mathematical sciences 15,1(2017),237–259.

技术实现思路

1、本发明的目的在于基于现有技术的现状，提出一种基于gpu并行和路径拟合技术的路径积分随机方法，用于gpu并行加速求解三维芯片中带混合边界条件的稳态热分析问题。

2、现有技术中，由于随机法求解结果一般具有的方差，其中n是路径数量[12]，因此需要超过数千条独立路径才能达到高精度。若能利用并行技术，则可大幅提升程序计算效率。然而，由于第三类robin边界在芯片热分析问题中普遍存在，通常每条随机行走路径需要执行数百万次跳跃才能确保所有路径最终温度收敛。这一特性对pirw程序的并行效果带来限制，即便所有路径全部并行，最终并行执行时间还是无法有效缩减。

3、本法利用gpu并行和拟合技术，加速基于路径积分随机方法(pathintegralrandom walk,pirw)求解三维芯片稳态热分析中的偏微分方程(partial differentialequations，pde)，即泊松方程(poissonequation)的混合边值问题(mixedboundaryvalueproblems,mixedbvps)。该方法能充分发挥gpu计算资源多的优势，有效地加速处理含第二类和第三类边界条件的问题，特别是针对第三类robin边界相关的feynman-kac泛函衰减较慢的难题，利用拟合技术进一步提升并行效率。

4、具体的，本发明中，

5、首先，将三维芯片的稳态热分析问题，转化为带混合边界条件的poisson方程边值求解问题；然后，基于feynman-kac公式，利用标准反射布朗运动，得到三维芯片稳态热分析问题的概率解；第三，将所有输入参数传输至gpu，根据所需的采样路径数n启动n个gpu线程，并初始化程序；第四，依据当前点xt所在的位置，在gpu中执行相应的随机行走过程；在无源区和热源区分别使用wos和wog进行随机行走；针对三类不同的边界条件，分别采用直接终止、计算局部时间、计算局部时间并更新三种不同的计算策略；第五，判断当前路径的跳转步数是否超过间距跳转步数δ的整数k倍，即np≥kδ，若是，则转至第六步，否则转至第四步；第六，记录第i条路径按照feynman-kac公式的离散形式计算跳转步数为kδ时的温度累计值ti，kδ(x)，若所有n条路径都至少完成kδ步跳转，转至第七步，否则转至第四步；第七，计算n条路径跳转步数为kδ时第三类robin边界和热源区对累计温度值的总贡献的平均值加入此平均值后在温度拟合窗口w内拟合温度判断该拟合温度是否已经收敛，即判断若是，则终止程序并计算待求点x处的温度t(x)作为输出，否则k增加1，转至第四步。

6、更具体的，

7、本发明提出一种基于gpu并行和温度拟合的路径积分随机加速方法，用于快速求解三维芯片中带混合边界条件的稳态热分析问题，其流程图如图1所示，其包括：

8、输入参数：区域ω及其边界的几何信息，区域边界上三种边界γi和边界条件i＝1，2，3，第一类吸收边界附近带状区域的厚度第二类和第三类反射边界附近带状区域的厚度及带状区域内wos跳转半径δx，待求点x的位置，采样路径数n，启动拟合的间距跳转步数δ，判断拟合值的收敛标准εth，温度拟合窗口w；

9、输出结果：待求点x处的温度t(x)。

10、其中具体步骤包括：

11、步骤1：将三维芯片的稳态热分析问题，转化为带混合边界条件的poisson方程边值求解问题；

12、步骤2：基于feynman-kac公式，利用标准反射布朗运动，得到三维芯片稳态热分析问题的概率解；

13、步骤3：将所有输入参数传输至gpu，根据所需的采样路径数n启动n个gpu线程，并初始化程序；

14、步骤4：依据当前点xt所在的位置，在gpu中执行相应的随机行走过程；在无源区和热源区分别使用wos和wog进行随机行走；针对三类不同的边界条件，分别采用直接终止、计算局部时间、计算局部时间并更新三种不同的计算策略；

15、步骤5：判断当前路径的跳转步数是否超过间距跳转步数δ的整数k倍，即np≥kδ，若是，则转至步骤6，否则转至步骤4；

16、步骤6：记录第i条路径按照feynman-kac公式的离散形式计算跳转步数为kδ时的温度累计值ti,kδ(x)，若所有n条路径都至少完成kδ步跳转，转至步骤7，否则转至步骤4；

17、步骤7：计算n条路径跳转步数为kδ时第三类robin边界和热源区对累计温度值的总贡献的平均值加入此平均值后在温度拟合窗口w内拟合温度判断该拟合温度是否已经收敛，即判断若是，则终止程序并计算待求点x处的温度t(x)作为输出，否则k增加1，转至步骤4。

18、本发明步骤1中，将三维芯片稳态热分析问题转化为带混合边界条件的poisson方程边值求解问题。

19、具体的，本发明方法解决的问题如图2所示。稳态热分析的三维区域为在区域ω内温度场t(x)满足poisson方程，区域边界为混合边界，即包含第一类dirichlet边界γ1、第二类neumann边界γ2和第三类robin边界γ3。给定三种类型边界条件，场域微分方程和边界条件分别为：

20、

21、其中，t(x)是区域ω内的温度，p(x)是热源的功率密度，k(x)为材料的导热系数。在第一类dirichlet边界γ1上，温度t(x)为已知的φ1(x)；在第二类neumann边界γ2上，温度t在边界处的外法向导数为已知的φ2(x)，在稳态热分析中，通常为绝热边界条件，即第三类robin边界γ3又称对流边界，对应于与周围流体交换热量的界面，例如冷却风扇和散热器的表面，参数其中h(x)为第三类robin边界处的热对流系数，φ3(x)为给定值。

22、本发明步骤2中，基于feynman-kac公式，利用标准反射布朗运动，得到三维芯片稳态热分析问题的概率解。

23、具体的，考虑一个满足(1)的区域，依据feynman-kac公式，该混合边值问题的概率解，即点x的温度值t(x)为[9]：

24、

25、其中，ex为数学期望，xt是反射布朗运动，l(t)是对应的局部时间(local time)，是反射布朗运动到达第一类dirichlet边界的时刻，是反射布朗运动到达第一类dirichlet边界的位置。(2)-(4)中的积分项分别表示第一类dirichlet边界条件，第二类neumann和第三类robin边界条件上，poisson方程右端项对于待测点温度的贡献。

26、feynman-kac泛函定义为[16]：

27、

28、反射布朗运动及对应的局部时间通过skorokhod问题基于布朗运动(brownianmotion，bm)定义得到。假设d是r3中一个具有c1边界的域，skorokhod问题的定义为：设f∈c([0，∞)，r3)，是从[0，∞)到r3的连续函数。一对(ξt，lt)是skorokhod方程s(f；d)的解，如果：

29、1.ξ在中是连续的；

30、2.l(t)是一个非减函数，当且仅当时增加，即

31、

32、3.skorokhod方程成立：

33、

34、其中n(x)表示处的单位外法向量,指标函数定义为：

35、

36、假设f(t)是从开始的布朗运动路径，并且(xt，lt)是skorokhod问题s(f；d)的解，则xt是从x开始的反射布朗运动。局部时间不是一个独立的过程，而是与标准反射布朗运动xt相关，定义为[17]：

37、

38、其中，d∈是包含的宽度为∈的带状区域并且l(t)是xt的局部时间，衡量反射布朗运动xt在[0,t]时间段内在第二类neumann或第三类robin反射边界的消失邻域d∈中耗费的时间量，这里

39、因此，feynman-kac公式提供了利用反射布朗运动(reflecting brownianmotion，rbm)xt求解poisson方程混合边值问题的方法，特别是包含第二类neumann边界或第三类robin边界的局部概率解，该概率解为对反射布朗运动xt进行路径积分的期望。其中，(2)-(4)分别是路径上关于第一类dirichlet边界条件、第二类neumann边界条件、第三类robin边界条件和泊松方程右端项的积分项。

40、通常，反射布朗运动采用wos方法进行模拟。反射布朗运动的起点x为球的中心，由球体的格林函数可知，反射布朗运动下一个点的位置是在球面上的均匀分布。因此，可以使用一系列球上的行走或跳转来模拟反射布朗运动的路径。同时，为了加快wos的过程，每一次跳转的球体半径被设置为最大值，即当前点到边界的距离。当反射布朗运动的点碰撞到第一类dirichlet边界时，该点会被边界吸收，反射布朗运动终止；当反射布朗运动的点碰撞到第二类neumann边界或第三类robin边界时，该点会被边界向内部反射，具体局部时间的计算细节参见步骤4.3。

41、本发明步骤3中，将所有输入参数传输至gpu，根据所需的采样路径数n启动n个gpu线程，并初始化程序。

42、具体的，gpu(graphics processingunit)是一种专为并行计算设计的处理器。以nvidia gpu为例，此类gpu采用支持cuda(compute unified devicearchitecture)编程模型的架构，其核心优势在于拥有数千个流处理器，能够同时执行大量线程任务。与cpu侧重复杂逻辑和低延迟不同，gpu采用simt(single instructionmultiplethreads，单指令多线程)执行模型，通过将计算任务拆分为大量可并行的子任务，实现数据级并行。例如在处理图像像素或矩阵运算时，每个流处理器可独立处理不同数据单元，显著提升吞吐量。

43、gpu并行计算通过层次化线程组织实现：线程组成线程块，多个线程块构成网格。硬件层面，流多处理器(streammultiprocessor,sm)动态调度线程束(warp)执行相同指令流，借助寄存器、共享内存等多级存储结构缓解访存瓶颈。这种架构特别适合处理具有高并行性、低数据依赖性的计算任务。

44、在启动gpu随机法求解程序时，指定总线程数为需要的路径数n，每条路径对应到gpu上的某一个线程。在第i条线程中，用若干变量保存该线程对应的路径的当前信息，包括当前跳转步数np、当前feynman-kac泛函的值、当前位置xt。第i条线程初始化程序时，将其当前跳转步数xt初始化为0，当前feynman-kac泛函的值初始化为1，当前所在位置x0初始化为x。另外，初始化程序时，所有路径已完成的δ整数倍跳转步数点k初始化为0，在gpu全局内存中预留空间用于保存跳转步数超过kδ时第三类robin边界和热源区贡献的累计温度值并将其中数据初始化为0。

45、本发明步骤4中，依据当前点xt所在的位置，在gpu中执行相应的随机行走过程；在无源区和热源区分别使用wos和wog进行随机行走；针对三类不同的边界条件，分别采用直接终止、计算局部时间、计算局部时间并更新三种不同的计算策略，并根据概率解的离散形式累计当前路径的温度值；

46、具体的，若当前点在无源区，则使用wos进行随机行走，即跳转步骤4.1；若当前点在第一类dirichlet边界区域内，则随机行走终止，即跳转步骤4.2；若当前点在第二类neumann边界区域内，则计算局部时间并累计温度值，即跳转步骤4.3；若当前点xt在第三类robin边界区域内，则计算局部时间，更新并累计温度值，即跳转步骤4.4；若当前点在热源区内，则采用wog进行随机行走并累计温度值，即跳转步骤4.5。步骤4.1到步骤4.5中所有跳转规则如图3所示。

47、由于较新的gpu架构(如volta架构)支持每个线程块中各线程的独立调度，所有的位置判断和温度累计都在同一个gpu函数中而不是被拆分为多个gpu函数模块[11]进行，以避免线程间通过gpu上的全局内存通信造成开销，如图4所示。

48、本发明步骤4.1若当前点xt在无源区，则使用wos进行随机行走。

49、具体的，若当前点xt在无源区且不在第一类dirichlet边界区域或第二类neumann边界、第三类robin边界的区域内，则使用wos进行随机行走，其跳转半径等于当前点xt到三类边界及热源区的最短距离。在以该距离为半径的球面上均匀采样，得到布朗运动的下一个随机点xt+1作为下一步的当前点xt，并使当前路径的跳转步数np增加1。

50、本发明步骤4.2若当前点xt在第一类dirichlet边界区域内，则随机行走终止。

51、具体的，在第一类dirichlet边界附近构建了一层极薄的带状区域其宽度为若当前点xt在第一类dirichlet边界区域内，则该点会被第一类dirichlet边界吸收到边界上，本条路径随机行走终止。

52、此终止时刻记为点xt距离第一类dirichlet边界最近的点记为该随机行走终止对温度贡献的积分项为：

53、

54、当不考虑第三类robin边界时，feynman-kac泛函可以忽略。然而，一旦引入第三类robin边界，此时需要在随机跳转路径中记录并更新的具体计算细节在步骤4.4中阐述。

55、需要指出的是，本发明提出的拟合技术加速不考虑终止于第一类dirichlet边界的路径。

56、本发明步骤4.3中若当前点xt在第二类neumann边界区域内，则计算局部时间并累计温度值。

57、具体的，在第二类neumann边界附近构建了一层极薄的带状区域其宽度为若当前点xt在第二类neumann边界区域内，需要计算每一点对应的局部时间。

58、若当前点xt在第二类neumann边界区域内，仍使用wos进行随机行走，跳转半径被固定为常数δx或2δx，每次跳转得到布朗运动的下一个随机点xt+1作为下一步的当前点xt，并使当前路径跳转次数np增加1。记当前点xt到第二类neumann边界的距离rr，当0＜rr＜δx时，令跳转半径r＝2δx，否则跳转半径r＝δx。采用这种固定步长的wos跳转方法，局部时间可表示为[18][19]：

59、

60、其中，δx是预设的在带状区域内的固定跳转步长，和分别为跳转点进出时的值，是路径从开始到时间tj停留在带状区域内的累积步数；值得注意的是，只有xt的跳转路径位于带状区域内才会对有贡献；

61、在区域中，wos的跳转步长为δx或2δx，因此跳转对应的局部时间为(δx)2/3或(2δx)2/3；因此，在实现中，本发明方法将视为一个递增的指标值；wos的每一步之后，的值相对于前一步的增量将为0、1或4，分别对应于随机点不在区域，在区域中并且跳转半径为δx，或在区域中并且跳转半径为2δx；其对温度贡献的积分项为：

62、

63、同样，当不考虑第三类robin边界时，泛函然而，一旦引入robin边界，此时，式(13)中的计算细节在步骤4.4中阐述。

64、本发明步骤4.4中，若当前点xt在第三类robin边界区域内，则计算局部时间，更新并累计温度值。

65、具体的，在第三类robin边界附近构建了一层极薄的带状区域若当前点xt在第三类robin边界区域内，则用wos进行随机行走，每次跳转得到布朗运动的下一个随机点xt+1作为下一步的当前点xt，并使当前路径跳转次数np增加1，具体方法与步骤4.3中的相同。

66、此外，还需要对进行更新。该更新与路径在区域内与第三类robin边界的碰撞情况有关。随机点在区域内随机行走时，经常未与边界碰撞就走出区域。本发明方法跟踪随机点进出区域的每个入口时刻tj，in和出口时刻tj，out，并计算在这个时间段内发生的每次碰撞对应的局部时间增量dl(t)。为了计算与碰撞点对应的需要记录wos路径进入和退出区域之间的系数增量以及同一时间段内路径与边界发生的碰撞次数mj，令

67、

68、然后，可以近似计算为：

69、

70、进而，其对温度贡献的积分项为：

71、

72、由于参数的值仅出现在第三类robin边界上，所以的值仅在时发生变化。当时，随机行走路径到达第一类dirichlet边界并被吸收，而表示随机行走结束时随机过程中第三类robin边界对路径的最终衰减。

73、feynman-kac泛函是第三类robin边界对路径造成的衰减，也就是说，对于只有第一类dirichlet边界和第二类neumann边界的混合边界问题，泛函等于1。但是，一旦涉及第三类robin边界，泛函将不会为1，并且会作为衰减系数影响所有积分项，不仅仅是与第三类robin边界条件相关的项。

74、因此，在对路径积分的计算过程中，函数需要作为一个变量动态更新。每次路径到达第三类robin边界时，的值会突然减小。并且，随着第三类robin边界反射的增加，的值逐渐从1变为0，导致后续积分项贡献的衰减和路径积分的最终收敛。值得注意的是，仅与第三类robin边界上的给定参数c相关，因此其值仅在路径到达第三类robin边界时才会发生变化。

75、本发明步骤4.5中，若当前点xt在热源区内，则使用wog进行随机行走，并累计温度值。

76、具体的，若当前点xt在热源区内，则使用wog进行随机行走，每次跳转得到下一个随机点xt+1作为下一步的当前点xt，使当前路径跳转次数np增加1，并累计温度值。为了实现wog，需对热源区域进行均匀网格离散，将热源区域划分为若干相同尺寸的长方体，得到热源区域的离散网格。其中，每个网格具有与相应的功率密度与功率值，网格节点即为网格的中心点。相邻两个网格节点之间可计算热阻和热导。如上下相邻的两个网格的长宽皆为hx和hy，高度分别为hzi和hzj，两个网格的热导率分别为ki和kj，则r(i，j)为两个网格节点之间的热阻，定义为：

77、

78、其中，g(i，j)为两个网格节点之间的热导。

79、从网格节点i跳转到临近网格节点j的概率p(i，j)为：

80、

81、其中degree(j)是网格节点i临近的网格节点数目。

82、xt在节点j获得的热源区温度贡献f(xt)在wog方法中可表示为r(j)，具体形式为：

83、

84、其中pj为在网格节点j处的功率。其对温度贡献的积分项为：

85、

86、在热源区与无源区交界处构建一层虚拟网格，其对应功率为零，进行wos和wog两种不同随机法的切换。当前点xt从无源区进入这层网格时，将会被拉到最近网格的中心点，然后采用wog方法随机跳转到一个临近网格节点或者无源区中。

87、本发明步骤5中，判断当前路径的跳转步数是否超过间距跳转步数δ的整数k倍，即np≥kδ，若是，则转至步骤6，否则转至步骤4；

88、具体的，第i条路径在每次完成步骤4中的跳转后，其对应的跳转次数np增加若干步。如果跳转次数np超过间距跳转步数δ的整数k倍，即np≥kδ，则应当记录这条路径当前的累计温度值，使转至步骤6，否则该条路径继续进行跳转，转至步骤4。

89、本发明步骤6中，记录第i条路径按照feynman-kac公式的离散形式计算跳转步数为kδ时的温度累计值ti，kδ(x)，若所有n条路径都至少完成kδ步跳转，转至步骤7，否则转至步骤4；

90、具体的，根据步骤4，按照feynman-kac公式的离散形式[13]，从点x出发的第i条路径在跳转步数为kδ时的实时累计温度值为：

91、

92、其中，其中j′表示路径中的每一个跳转点，表示在第二类neumann边界、第三类robin边界的区域中的跳转点，表示在热源区域内的跳转点。注意保存该温度值时，路径的实际跳转步数np可能略大于kδ，但计算时只计算从0步跳转到kδ步时的累计温度值。该温度值ti，kδ(x)将被保存在gpu全局内存中。完成保存后，判断是否所有路径均已完成kδ次跳转，若是，所有kδ步的累计温度采样已完成，可以计算平均值，转至步骤7，否则暂时不计算平均值，当前路径继续跳转，转至步骤4。

93、本发明步骤7中，计算n条路径跳转步数为kδ时累计温度值的第三类robin边界和热源区对累计温度值的总贡献的平均值加入此平均值后在温度拟合窗口w内拟合温度判断该拟合温度是否已经收敛，即判断若是，则终止程序并计算待求点x处的温度t(x)作为输出，否则k增加1，转至步骤4。

94、具体的，在实际三维芯片热分析中，通常存在大量的第二类neumann边界和第三类robin边界。对于这两类边界，由于feynman-kac泛函每条路径都需要大量的跳转步数才能使累计的温度收敛，这导致程序的运行时间较长。然而，在实际热分析应用中，当跳转步数足够大时，pirw获得的温度通常会开始呈指数收敛趋势，这一特性可以通过拟合来减少达到所需精度所需的步骤数。

95、由于问题中存在第一类dirichlet吸收边界时，路径通常快速终止，不会影响并行程序的执行时间，因此，本发明在拟合中不考虑第一类吸收边界。此外，在实际芯片热分析问题中，第二类neumann边界条件通常为绝热边界条件，其对累计温度值的贡献为0，不在拟合中考虑。因此，n条路径在kδ步时第三类robin边界和热源区对累计温度值的总贡献的平均值可以写为：

96、

97、其中i(.)为指示函数。

98、对公式(22)作近似，可以得到：

99、

100、

101、根据布朗运动的性质，在有限的空间内，经过足够多的跳转步数之后，所有路径的当前位置的空间分布将趋于稳定。由于只与当前步数j下所有路径当前位置的空间分布有关，经过足够多的跳转步数之后，会稳定为一个常数。又由于feynman-kac泛函的衰减由指数中的和dl(tj)控制，而和dl(tj)也只与当前步数j下所有路径当前位置的空间分布有关，因此，feynman-kac泛函会随着跳转步数的增加呈指数衰减。综合以上两点，当跳转步数足够多时，每一步的温度增量应该近似与当前的feynman-kac泛函成比例，即随着跳转步数的增加呈指数衰减，进而导致累计温度值开始呈指数收敛，如图5所示。

102、基于上述假设，在n条路径都已经完成kδ步跳转时，使用在温度拟合窗口w内多个跳转步数点处第三类robin边界和热源区对累计温度值的总贡献的平均值通过公式(25)拟合估计的收敛温度tc和参数a、b和c。拟合得到的tc记为

103、

104、实际运行中，通过判断拟合值是否收敛来决定是否结束程序。将当前拟合温度与前一个拟合温度进行比较，若则用最新的拟合值作为最终结果并终止所有路径。

105、若拟合值已经稳定，即则终止程序，并输出点x的温度

106、

107、其中，为终止时刻，为终止时点xt距离第一类dirichlet边界最近的点，j′表示路径中的每一个跳转点，表示在第二类反射边界区域内的跳转点。

108、否则，即k增加1，转至步骤4继续进行随机行走。

109、本发明方法的优点在于：

110、(1)在gpu上实现并行高精度随机法三维芯片稳态热分析方法，利用gpu上数千个计算核心并行计算，显著减少了pirw程序的计算时间。

111、(2)提出一种拟合技术，利用温度值近似以指数方式收敛的特性，减少随机法中路径运行至温度收敛所需的步数。

112、本技术针对三维芯片中混合边界条件下稳态热分析问题，提出的gpu并行加速方法是求解laplace方程和poisson方程混合边值问题的有力工具，对于解决大规模三维芯片复杂热问题具有显著的速度优势。并且，本方法可以拓展到其他工程领域中具有混合边界条件的laplace方程和poisson方程的场分析问题中。

113、本方法的优点有：

114、(1)在gpu上实现并行、高精度随机法三维芯片稳态热分析方法，利用gpu上数千个计算核心并行计算，显著减少pirw程序计算时间。

115、(2)利用随机行走中温度值近似以指数方式收敛的特性，提出一种拟合技术实现随机行走提前终止，减少随机法中每条路径运行至温度收敛所需的步数，以达到加速的目的。