PODTherm-GP: A Physics-based Data-Driven Approach for Effective Architecture-Level Thermal Simulation of Multi-Core CPUs

To obtain detailed thermal profiles in semiconductor chips, usually the FEM, FVM or FDM is applied. For example, 3D thermal ADI based on the FDM was developed by Wang et al. [16]. Xu [17] used ANSYS based on the FEM to perform thermal simulation for a quad-core chip. COMSOL was used by Vaddina et al.  [18] to conduct thermal modeling and analysis of 3D stacked structures. With rapid advance in multi-core processors [19], localized power density has been significantly enhanced and it is desirable to attain high enough resolution in a thermal prediction to capture small-size hot spots. Although, DNS offers high thermal resolution and high accuracy, it is impractical to apply DNS at the chip architecture level due to their computationally intensive nature.

Lumped element thermal circuits have been widely used at the architecture level due to their simplicity and efficiency [20]. For example, HotSpot [21] and 3D-ICE [13] are some of the popular thermal-circuit architecture level simulators. HotSpot has been applied by Glocker et al.  [22] to study the temperature distribution in a 16-core system and by Florea et al. [23] to provide thermal analysis for enhancing the Sniper multicore simulator. HotSpot and 3D-ICE have been integrated into the performance-power-thermal simulation toolchains, such as HotSniper  [24], HotGauge  [25] and CoMeT  [26]. The efficiency of the thermal circuit model is however accomplished by using large lumped elements that are not able to adequately account for distributed heat transfer and may lead to an inaccurate prediction [27]. To improve the accuracy, the size of the RC elements needs to be reduced substantially, as done in the grid model of HotSpot [28] that is then similar to the computationally intensive FDM.

In the conventional Green’s function approach, the single Green’s function is calculated in response to a unit point heat source at the center of a large chip, assuming no boundary influence. The Green’s function is a spatial impulse response of the chip, and the thermal solution is constructed by superposition of the impulse responses to the point sources at different locations. It is thus inherently difficult to include boundary conditions (BCs) in thermal simulation of a finite domain [29, 30] or to use Green’s function in situations where the power source is close to the edge of the chip [31]. It is also difficult to apply the conventional approach to transient thermal simulation [29, 32]. The method is however significantly more efficient than the FEM, FVM or FDM because it only considers a single layer where the power sources are generated. Various efforts have been made to overcome the aforementioned limitations by implementing different techniques with some significant revisions, which inevitably reduce computing efficiency. For example, the Green’s function is applied to correct the corner and edge effects based on the method of image for the adiabatic BC [33]. This is however not valid for different types of BCs. Using the power blurring method, together with the spatiotemporal pulse, transient thermal simulation becomes possible within the Green’s function framework [33]. These however still offer 2D temperature on the power dissipation layer. An effort has been made to include 3D temperature profile with multilayer Green’s functions [34]; it is however limited to steady state simulation.

Machine learning methods have recently been applied to investigate various aspects in semiconductor chips, including multi-core chips, [35, 36, 37]. Among different machine learning techniques, the neural network methods have been the most common approaches for thermal simulation of semiconductor chips due to their simplicity and efficiency [35, 36]. Using neural networks, thermal data are needed to train each thermal cell (or neuron) in the domain to respond to heat excitations (heat sources and BCs) and neighboring neurons. Memory space and computational time needed in thermal simulation are thus determined by the selected number of neurons and the selected distance for close connections of the individual neurons. Although these approaches are considerably more efficient than the DNSs; they are in general limited to the thermal solution over a small number of nodes, instead of the detailed thermal profile of the chip. In addition, the approaches are based on data statistics without being guided by fundamental physical principles and thus in general leads to erroneous solution in case of extrapolation.

A different learning algorithm derived from a projection-based reduced order model has been investigated in the past based on POD [38, 39]. A tutorial for POD is included in  [40] that offers an intuitive and easy-to-follow presentation about the basic concepts and approach. The POD projects the thermal problem from a physical domain to a functional space to reduce the DoF. Instead of assuming its basis functions as done in many other projection-based approaches, such as the Fourier series, Legendre polynomial, Bessel functions, etc., the basis functions (or POD modes) that constitute the POD space are generated (or trained) via thermal solution data for the problem of interest. The thermal data are collected from accurate DNSs for the trained POD modes to capture the essential thermal information induced by the parametric variations, including heating power and BCs. The POD offers an optimal set of basis functions that lead to a smallest least square (LS) error limited by the quality of the thermal data. To incorporate the physics principles of heat transfer in POD simulation, the dynamic heat transfer equation for the simulation domain is further projected onto the POD space. This rigorous procedure guides the POD thermal solution to comply with the dynamic heat transfer equation. Cheng’s group has applied the POD simulation methodology to develop efficient and accurate thermal simulation models for SOI and FinFET devices, FinFET circuits at the gate level and interconnects [8].

The major advantage of the POD simulation methodology, guided by fundamental physics principle, is to reduce a large-dimension problem to an extremely small dimension with only a handful of modes while maintaining the accuracy comparable to and the resolution equivalent to DNS. A post process is however needed to reconstruct the dynamic temperature in the physical domain, which is the major computational bottleneck. Unlike the DNS where the entire simulation time and domain needs to be included in the simulation, the post processing calculations in the POD approach for realistic situations only need to be carried out in small high-temperature areas of the chip at each selected interval of time. This will significantly minimize the computing time and memory space needed in the post process.

The POD modes are optimized by maximizing the mean square inner product of the thermal solution data with the modes over the entire domain [38, 39], subjected to dynamic or static parametric variations. In our study, the thermal data in space is collected at each simulation time step in the DNS. This maximization process leads to an eigenvalue problem described by the Fredholm equation,

where $\lambda$ is the eigenvalue corresponding to the POD mode (i.e., eigenfunction $\varphi$) and $R(\vec{r},\vec{r}^{\prime})$ is a two-point correlation tensor given as

with $\otimes$ as the tensor operator and the angle brackets $\langle\rangle$ indicate the average over a number of thermal data sets. After solving the eigenvalue problem in (1), temperature can then be represented by

where $a_{i}$ are the weighting coefficients for $\varphi_{i}$ and $M$ is the number of POD modes selected to reconstruct the temperature solution. If the eigenvalue of the data decreases rapidly for the higher modes, only a small DoF (or $M$) is needed to reach an accurate thermal prediction.

In order to predict the dynamic thermal distribution given in (3) in a domain structure, the coefficients $a_{i}(t)$ need to be determined. To achieve this, the heat transfer equation is projected onto the POD space using the Galerkin projection,

where $k$, $\rho$ and $C$ are the thermal conductivity, density and specific heat, respectively, $P_{d}(\vec{r},t)$ is the interior power density, $S$ is the boundary surface and $\vec{n}$ is the outward normal vector on the boundary surface. With the selected POD modes,  (4) can be rewritten as an $M$-dimensional ordinary differential equations (ODEs) for $a_{i}(t)$,

where $c_{i,j}$ are the elements of the thermal capacitance matrix in the POD space and they are defined as

$g_{i,j}$ and $p_{j}$ are the elements of the thermal conductance matrix and power vector in the POD space, respectively. In  (5), $g_{i,j}$ and $p_{j}$ may also vary with BCs. In this work, two kinds of BCs, adiabatic and convective BCs, are applied to the boundary surfaces of the chip. For the adiabatic surface, the heat flux is zero, and $g_{i,j}$ and $p_{j}$ can be written as

For the convective boundary, the heat flux on the surface is described by

where $h$ is the heat transfer coefficient and $T_{amb}$ is ambient temperature. Using (8) in (4), $g_{i,j}$ for the convective BC are given by

and $p_{j}$ becomes

Once the dynamic power consumption is obtained from gem5 and McPAT (described in Sec. IV), the interior power source strength in POD space given in (7) and (10) can be pre-evaluated. The BC of the substrate bottom is modeled by convective heat transfer given in (8) with a heat transfer coefficient determined by structure dimensions and material properties of heat spreader and heat sink with $T_{amb}$ = 45 °C. All other boundaries are adiabatic. The coefficients in (5) are all pre-evaluated once the modes are determined. With $a_{i}(t)$ solved from (5) in POD simulation, the temperature solution can be predicted from (3).

The simulation domain of the selected multi-core CPU, whose floorplan is shown in Fig. 2 [41], includes a top heating layer with a thickness of $55.8~{}um$ and a substrate with a thickness of $241.8~{}um$. The heating layer covers the top thickness where power is dissipated (such as the layers of devices and interconnects). Convective heat transfer described by a heat transfer coefficient given in (8) is applied to the bottom of the substrate. The heat transfer coefficient is modeled by thermal resistances of the heat spreader, thermal interface material and heat sink calculated from the chip dimensions and material properties, similar to other studies [21, 28]. The remaining boundary surfaces are assumed adiabatic. The dynamic power in each functional unit generated from gem5 and McPAT is averaged over 10k CPU cycles at 2.3 GHz. To investigate the effect of the training data quality on robustness of the POD model, two sets of thermal data are collected from FEniCS-FEM with meshes of $256\times 256\times 14$ and $1104\times 620\times 14$, or spatial resolutions of $0.121\times 0.083\times 0.023~{}mm^{3}$ and $0.028\times 0.035\times 0.023~{}mm^{3}$ in $x$, $y$ and $z$, respectively, using the same dynamic power map. The POD models built upon these 2 data sets with coarser and finer resolutions are referred to as POD Model-A and Model-B, respectively. It is noted that the coarser mesh used in this study is finer than those used in most of studies for chip-level thermal simulations, unless an extremely fine resolution is needed  [33, 31].

The POD modes and eigenvalues are determined from (1) via the method of snapshots [45, 8] using the training data collected from FEniCS-FEM simulations. The POD modes are thus tailored to essential information embedded in the training data to account for variations of BCs and the dynamic power map. Each eigenvalue represents the mean squared temperature variations captured by its mode and therefore reveals the importance of the mode. The number of POD modes needed in (3) for an accurate prediction can then be estimated from

where $Err_{Theo}$ is the theoretical LS error and $N_{s}$ is the number of snapshots or sampled data sets. (11) offers an ideal LS error but the LS error resulting from PODTherm-GP is usually larger than $Err_{Theo}$ due to the limitations of numerical accuracy and computer precision.

Fig. 3 describes the eigenvalues of the collected data for Model-A and Model-B in descending order. Since the resolution of the collected data sets are fine enough, the eigenvalues for these 2 POD models are nearly identical. The eigenvalue decreases significantly in the first several modes and shows a slower decreasing rate in the higher modes. It eventually becomes nearly flattened beyond 170 modes after a reduction by 16 orders of magnitude due to the limit of computer precision. The zoom-in spectrum in the inset shows that the eigenvalue drops more than two orders of magnitude from the first to the second mode and nearly four orders to the third mode. Based on  (11), it is expected that either Model-A or Model-B is able to offer an accurate thermal prediction with 3 or more modes if the quality of the training data collected by the FEniCS-FEM is adequate.

To construct a POD model, its model parameters in the projected ODEs in (5) need to be pre-calculated using its POD modes $\varphi_{i}(\vec{r})$ of the training data, as defined in  (6)-(10). PODTherm-GP thus solves the ODEs for $\vec{a}(t)$ with a selected number of modes $M$, and a post process via  (3) is then needed to calculate the dynamic thermal distribution in the CPU. Using the POD modes with two different resolutions, the following demonstrations focus on validating the accuracy and efficiency of the POD simulation technique and examining how the resolution impacts the quality of the training data and the generated modes and thus the model accuracy. A validation is also performed beyond the training range to examine the robustness of PODTherm-GP in the case of extrapolation.

In the demonstrations, each of Model-A and Model-B is applied to thermal simulation of the selected multi-core CPU, compared to the FEniCS-FEM simulation with the spatial resolution identical to that of its POD modes. The training is carried out over a simulation time of 4.35 $ms$ for both Model-A and Model-B, and the POD models are validated against FEniCS-FEM for 4.35 $ms$ and 6.09 $ms$. In any comparison between simulations of a POD model and FEniCS-FEM, BCs and dynamic power map are identical. However, the dynamic power map used for the demonstration is different from that used for the collection of training data, as discussed in Sec. IV.

Fig. 4 illustrates the LS errors of both POD Model-A and Model-B for the thermal simulations of the selected multi-core CPU, with respect to FEniCS-FEM simulations. The numerical LS error of the POD model is estimated over the entire simulation time from

where $T_{i}(\vec{r})$ and $e_{i}(\vec{r})$ are temperature given by FEniCS-FEM and the temperature difference between the FEniCS-FEM and the POD model at $i$-th time step, respectively, and $N_{t}$ is the total number of time steps.

As can be seen, the theoretical LS errors for Model-A and Model-B basically overlap due to the nearly identical eigenvalue spectrum shown in Fig. 3. The theoretical LS error decreases rapidly from 5.1% for one mode to 0.71% with 3 modes. When using 3 modes, the numerical LS error shown in Fig. 4(a) in the entire chip however reach 2.9% for Model-A (coarser resolution) and 2.1% Model-B (finer resolution) for the 4.35 $ms$ simulation. Beyond 3 modes, the error from the coarser mesh model remains nearly unchanged but from the finer mesh model continues decreasing slowly to 1.93%. For the heating layer in Fig. 4(b), the LS error for the 4.35 $ms$ simulation with 3 modes using Model-A or Model-B becomes 1.1% or 0.76%, respectively. The error from Model-A continues decreasing to 0.91% with 8 modes while the error from Model-B decreases to 0.49% with 8 modes and stays near 0.51% beyond 8 modes. The errors in the entire chip and the heating layer are clearly improved when the POD modes are generated from the finer mesh training data.

For the 6.09 $ms$ case with simulation time beyond the training range, as expected, the LS error of the POD model increases. However, Figs. 4(a) and (b) clearly show that, when more modes are included in this extrapolation case, the LS error of the 6.09 $ms$ simulation declining toward the error of the 4.35 $ms$ case. For example, when using 3-11 modes, the LS error of Model-B for the entire chip remains near 1.93%-2.1% for the 4.35 $ms$ case, and its error for the 6.09 $ms$ case is equal to 2.5% with 3 modes and reduces to 2.2% and 2.09% with 7 and 11 modes, respectively. It is observed as well for Model-A in Fig. 4(a) that the 6.09 $ms$ curve declines with more modes and these 2 LS-error curves nearly merge with 10 or more modes. The similar phenomenon is also shown in Fig. 4(b) for the heating layer, when comparing the LS errors induced by these two POD models using more modes.

As demonstrated in the 4.35 $ms$ case, for Model-A with 3 modes compare to FEniCS-FEM, a reduction in the DoF over 5 orders of magnitude ($256\times 256\times 14$ vs. $3$) with a LS error of 2.9% for the entire chip or 1.1% in the heating layer can be achieved. For Model-B with 3 modes, a reduction of DoF more than 6 orders of magnitude ($1104\times 620\times 14$ vs. $3$) with an LS error of 2.1% or 0.76% for the entire chip or the heating layer, respective, can be reached. In this 4.35 $ms$ case, if 7 modes are used, an LS error as small as 1.95% or 0.52% can be accomplished over the entire chip or in the heating layer, respectively, and the reduction in the DoF is still as large as 6 orders.

To understand the insight into the LS errors presented in Fig. 4, how the dynamic thermal profile influenced by the accuracy of the POD modes resulting from the quality of the training data is examined below in detail.

The temperature evolution at the intersection of Paths A and B predicted by POD Model-A and Model-B is compared against the FEniCS-FEM result in Figs. 5(a) and (b), respectively, over a simulation time of 6.09 $ms$. At this particular location, both POD-Galerkin models with 3 or more modes offer a very good agreement with FEniCS-FEM even beyond the training time (4.35 $ms$). When using 7 modes, the absolute percent errors for $t>2.2~{}ms$ become nearly flattened and stay near or below 1.62% from Model-A and 0.85% from Model-B, where the finer-resolution model leads to a higher accuracy. Because $T_{amb}$ is used as the initial chip temperature, the percentage error w.r.t. $T_{amb}$ looks evidently large for $t<0.55~{}ms$ even though the absolute error is less than $0.1^{\circ}$C at $t\approx 0.55~{}ms$ and much less than $0.01^{\circ}$C near $t=0^{+}~{}ms$.

The temperature profiles along Paths A and B are shown in Fig. 6 for Model-A and in Fig. 7 for Model B. Based on the CPU floorplan in Fig. 2, Figs. 6 and  7 show that temperatures in Cores 12, 16 and 17 appear to be higher than other functional units of the multi-core CPU resulting from the core assignment based on the selected benchmarks and their computational intensities. In general, POD Model-A with 3 modes offers an accurate thermal prediction along Paths A and B in the chip, as presented in Fig. 6. However, in the high temperature regions, Cores 16 and 17, there is a 2%-3% deviation from the FEniCS-FEM result. When a finer mesh is used in FEniCS-FEM to collect the training data to generate the POD modes, Model-B with 3 or more modes offers an excellent thermal prediction along both Paths A and B, as shown in Fig. 7. As also illustrated, the error derived from Model-B (a finer mesh in Fig. 7) is 50% less than that from Model-A (a coarser mesh in Fig. 6) in nearly all locations along Paths A and B.

To understand how the mesh size impacts the POD model accuracy, the profiles of POD modes and gradients of the modes for Model-A and Model-B are examined in Fig. 8. The difference of the mode gradients between Model-A and Model-B along Paths A and B are also included. Because the eigenvalue spectrums for these 2 POD models are nearly identical, as given in Fig. 3, their modes are very close. However, a close look at the slopes of the modes in each direction reveals the hidden problem of the model built upon the coarser mesh. As shown in Fig. 8, the gradient differences of the POD modes between Model-A and Model-B at some locations with larger thermal gradients is relatively large. For example, as shown in Figs. 8(a)-8(c), a slope difference of the POD modes in $x$ around 25%, 75% or 45% is observed in Modes 1, 2 or 3, respectively, near $x=4.9~{}mm$ where the highest thermal gradient in $x$ is observed in Figs. 6(a) and  7(a). Evident slope differences of the modes in $x$ are also observed in other large slope locations near $x=12.4~{}mm$ and $16.15~{}mm$. Along Path B, the slope difference in $y$ shown in Figs. 8(d)-8(f) is smaller but still evident for the first 2 or 3 modes at $y=10.4~{}mm$ and $17.9~{}mm$ but the difference at $13.4~{}mm$ is observed only in the third mode.

Accuracy of the POD model in (5) is dictated by the quality of the coefficients evaluated in (6)-(10). The thermal conductance elements $g_{i,j}$ in the POD space results from the projection of the heat flux ($-k\nabla T$) onto the POD modes and thus strongly depends on gradients of the modes, as defined in (7) and  (9). The accuracy of the POD model is thus significantly influenced by the gradients of the generated POD modes that are directly impacted by the quality of the training data. This explains why Model-B built upon the thermal data collected from the FEniCS-FEM simulation with a finer mesh gives rise to a considerably smaller LS error than Model-A, as discussed above in Fig. 4.

To further examine the robustness of PODTherm-GP in the extrapolation case, the thermal profiles at 6.09 $ms$ predicted by Model-A and Model-B are illustrated in Figs. 9 and 10, respectively, compared to FEniCS-FEM results. For the case with a coarser mesh in Fig. 9, the thermal profiles along Paths A and B from POD Model-A with 3 or 5 modes agree reasonably well with those from FEniCS-FEM, except for the location near $x<3~{}mm,19~{}mm<x<23~{}mm$, and $x>27~{}mm$ in Fig. 9(a) and $5~{}mm<y<10~{}mm$ in Fig. 9(b). Fig. 4(b) illustrates that, with 7 or more modes, Model-A starts improving its LS error for the simulation beyond the training period and a better agreement with FEniCS-FEM is observed in Figs. 9(a) and (b). Because of considerably higher power applied to Cores 3 and 8 and a longer simulation time than the training conditions, it requires 19 modes to significantly improve the prediction from Model-A. The absolute error included in Fig. 9 also reveals a significant improvement when using 19 modes, compared to 7 modes although an error as small as 1.5%-2% still remains near $13~{}mm<x<15~{}mm$ in Core 17 (similar to Fig. 6 (b)). For the finer mesh case in Fig. 10, a similar improvement is obtained for the Model-B simulation beyond the training period when more modes are used, and a considerably better prediction is observed for $x<3~{}mm$ and $13~{}mm<y<15~{}mm$, where an excellent agreement with FEniCS-FEM is observed.

Results shown in Figs. 4, 9 and  10 demonstrate that, even with core power and simulation time beyond the training conditions, PODTherm-GP offers a prediction with a reasonably good accuracy with just 7 modes and more modes are needed to reach a very accurate prediction. Such an extrapolation capability stems from the Galerkin projection of the heat transfer equation. Unlike most machine learning methods that usually lead to an unpredictable solution when simulation settings are outside the bounds of the training conditions, the Galerkin projection enforces the prediction to comply with the physical principles embedded in  (5). These results also confirm the significant influence of data quality (resulting from the mesh size in this case) on the prediction accuracy of PODTherm-GP in the extrapolation case.