I. INTRODUCTION

Artificial intelligence (AI) and machine learning (ML) have been widely used in most of domains of technology [1], [2]. The models used in AI and ML are trained by processing millions of crawled data giving rise to considerable demand for highperformance processors [1] To satisfy the need, more cores are integrated on a semiconductor chip, and the density of transistors and power dissipation have been increasing dramatically in recent years, which has led to high temperature and hot-spot generation due to severe joule heating. High temperature and hot spots contribute to not only degradation of performance but also deterioration of reliability [3], [4]. To reduce temperature and suppress hot spots in high-performance processors, the general practice is to apply effective thermalaware task scheduling and thermal management, which however requires effective and accurate chip-level thermal-simulation techniques.

Several approaches have been developed for the thermal simulation of semiconductor chips; each of them offers a different level of efficiency and accuracy. Among these approaches, direct numerical simulations (DNSs) based on either the finite element method (FEM) or finite difference method (FDM) provide accurate and detailed thermal analysis at the expense of a large number of degrees of freedom (DoF). Many open-source or commercial DNS tools are available for such applications, for example, FEniCS [5], ANSYS [6], COMSOL [7], etc. These DNSs, although offering accurate thermal solution with fine resolution, demand extensive computational resources and are impractical for chip-level thermal simulations.

To conduct the chip-level thermal simulation efficiently, the lumped RC thermal circuit model has been used to predict the thermal profile in large-scale semiconductor chips; for example, the block model of HotSpot [8]- [10] is one of most popular thermal simulators using the compact RC thermal model for chip-level thermal simulations. Due to the large RC lumped element, the RC thermal circuit model is not able to capture the small-size hot spots in semiconductor chips but only offers average temperatures for the large RC elements. With the approximation associated with large lumped element, heat flux at the element interfaces cannot be estimated accurately. The accuracy of the block model of HotSpot has been challenged due to the inaccurate thermal prediction for some floorplans, compared to DNS [11]. To improve the accuracy of the block model of HotSpot, the grid model of HotSpot [12] was developed, where smaller elements are allowed to provide a more detailed/accurate temperature prediction. However, when using very small elements for better accuracy, the grid model of HotSpot is equivalent to the FDM and becomes prohibitive for chip-level simulation.

To enhance the efficiency of chip-level thermal simulations, another strategy is to develop a spatial impulse response (or the Green's function) [13]- [15] of the selected chip. The Green's function is usually pre-trained by the thermal solution derived from DNS in response to a unit point heat source at the center of the chip. The spatial temperature solution is then obtained by a convolution of the pre-trained Green's function with the power profile. However, for the Green's function method, it is difficult to apply boundary conditions (BCs) [13], [14] or to perform transient thermal simulation [13], [15]. In addition, the training of the Green's function using DNS of the entire chip is extremely time consuming [13], especially if a high resolution is needed to capture the localized hot spots. As the technology node is further reduced and more cores are integrated on a chip, the computation of the Green's function is becoming more intensive and impractical for developing such a thermal model for the entire chip, especially when a high resolution is needed.

An alternative is to use a reduced-order simulation model enabled by a data-driven approach based on proper orthogonal decomposition (POD) [16], [17]. This approach projects a dynamic thermal problem from the physical domain onto a functional space (POD space) described by a finite set of basis functions (also called POD modes). To derive an optimal set of modes, dynamic thermal data accounting for parametric variations of interest, such as variations of heat excitations and BCs, are obtained from DNSs to train the POD modes. The POD model constructed by these trained robust modes is therefore able to respond accurately to the parametric variations within or near the training conditions with a very small number of DoF. In addition to the high accuracy and efficiency, the POD model also offers the temperature profile as detailed as DNS.

The POD simulation approach has been shown to be effective in many areas of research [18]- [27]including thermal simulations of integrated circuits and CPUs [20]- [22], [27]. However, similar to the problem encountered in the pre-training of the Green's function, a long simulation time and massive thermal data needed to train the POD modes become prohibitive for larger chips with high resolutions. To overcome the difficulty, the multi-block POD methodology is proposed for large-scale chips, such as multicore/manycore processors. In the multi-block POD model, the domain decomposition technique is implemented to partition a large semiconductor chip into smaller building blocks, such as cores, caches, I/O units, etc. For each small block, a set of POD modes and the model parameters can be generated more efficiently and stored into a technology library. The POD model for the entire chip can then be constructed by gluing these POD blocks with the discontinuous Galerkin (DG) method [28], [29]. This method is applied to stabilize the numerical solution at the interface by enforcing the heat flux continuity but allowing a small temperature discontinuity (i.e., the weak boundary condition) in an average sense at the interface between any 2 neighboring blocks. With the multi-block POD model for a large chip partitioned into a large number of building blocks, parallel computing can also be implemented in POD mode generation and thermal simulation to further enhance the computational efficiency.

Continuing a previous study [27], this work investigates a two-block POD model that projects two building blocks (Core 1 and its adjacent L2 cache) in AMD Athlon II X4 610e [30] to a POD space described by the 2-block POD modes. For each building block, DNSs are performed in FEniCS [5], an opensource FEM platform, to collect temperature data for the extraction of POD modes. The two-block POD model is demonstrated and verified against the DNS, and it has shown that the POD results are in very good agreement with the DNS with almost 4 orders reduction in the DoF.

II. THERMAL SIMULATION METHODOLOGY BASED ON POD

A. Single-block model

Using the POD method, the physical domain is projected onto a mathematical space represented by a finite number of POD modes. Temperature in space and time 𝑇(𝑟 ⃑, 𝑡) can then be represented by a linear combination of the selected POD modes 𝜑 ! as

where 𝜑 ! is the i-th POD mode, M is the number of selected POD modes which determines the accuracy and efficiency of the POD approach and 𝑎 ! (𝑡) is the time-dependent coefficient of the i-th POD mode.

To obtain an optimal set of the POD modes, each POD mode is obtained by maximizing the mean square inner product of the thermal solution with the modes via the following equation

where Ω is the physical domain of the selected structure and the brackets 〈 〉 denote the average over the collected thermal solution data. For dynamic thermal simulation, the average is computed over temporal samples (snapshots) obtained from DNSs. The maximization process in (2) gives rise to a Fredholm equation shown below for the POD modes,

where 𝐑(𝑟 ⃑, 𝑟 ⃑ ( ) is a two-point correlation tensor expressed as

With the temperature data 𝑇(𝑟 ⃑, 𝑡) of the simulation domain collected from DNSs, the method of snapshots [25], [26] is applied to solve the eigenvalue problem in (3) for the eigenvalues 𝜆 ! and POD modes 𝜑 ! .

With the generated POD modes, the heat conduction equation can be projected onto a POD space represented by the POD modes using the Galerkin projection,

surface of the selected domain and 𝑛 8⃑ is the outward normal vector of boundary surface. Substituting (1) into (5), it leads to an M-dimensional ordinary differential equation (ODE) for 𝑎 ! (𝑡),

where 𝑃 ! representing the last 2 terms of ( 5) for the i-th mode is the power density dissipated in the POD space and can be preevaluated since the shape of power density is predefined, and 𝑐 !,. and 𝑔 !,. are the elements of thermal capacitance and thermal conductance matrices in the POD space and defined as As presented above, the POD model development consists of thermal data collection from DNS, calculations of POD modes and eigenvalues from (3) using the snapshot method, and evaluations of model parameters in (7). This training process could be computationally intensive for a large simulation domain with a high resolution. To minimize the computational resources in the training, the large domain is partitioned into smaller building blocks, which is presented next.

B. Multi-block model

When placing block together, the last term of (5) needs to be reformulated to account for heat flux across the interface between adjacent blocks. The DG method [28], [29] is applied to properly enforce the interface thermal continuity, and (5) for the multi-block model becomes

where ⟦ * ⟧ and { * } indicate difference and average across interface, respectively, and 𝜇 is the penalty constant defined as 𝑁 0 𝑑𝑟 ⁄ (𝑑𝑟 is the size of the local element with Nµ as the penalty number). 𝑆 is the interface surface between two adjacent blocks. Nµ can be adjusted to balance discontinuities between temperature and heat flux at the interface to minimize the least square (LS) error and to stabilize the numerical solution.

For a two-block POD model including the heat flux exchanges via the interface, the matrix equation for both POD blocks becomes

where 𝐵 2,3 indicates the interface between 𝑝-th and 𝑞-th blocks. 𝑪 2 and 𝑮 2 are the thermal capacitance and thermal conductance matrices of 𝑝-th blocks and their elements are given by (7).

Compared to the thermal conductance matrices of single-block model, an extra thermal conductance matrix, 𝑮 !,# !,# , is included for 𝑝-th block to consider the effect of 𝑞-th block with respect to temperature and heat flux at the interface 𝐵 2,3 and is given as

where 𝑀 2 is the number of selected POD mode of 𝑝 -th block, and the element of 𝑮 !,# !,# is given by

In the matrix equation ( 9), the 𝑝-th block is coupled with its adjacent 𝑞 -th block via 𝑮 2,3 . If 𝑝 -th block is not directly adjacent with 𝑞-th block, 𝑮 2,3 = 𝟎. Otherwise, it is given as

where 𝑀 3 is the number of selected modes of 𝑞-th block and 𝑔 !,$,%,& is given by

III. CHIP-LEVEL THERMAL SIMULATION USING MULTI-BLOCK

POD MODEL The AMD ATHLON II X4 610e processor is selected in this investigation, which consists of four cores, two L2 caches, one northbridge, three I/O units and one DDR3 module, as shown in Fig. 1. The dimension of the quad-core chip is 14mm × 12mm × 242µm (length×width×thickness) and the material property is listed in Table I. In the single-block model, the thermal simulation is performed via DNSs over the entire chip to collect temperature data and it is, as discussed above, computationally intensive for a large chip with a high resolution. In the multi-block model, the temperature data is, however, independently collected for each building block in the entire domain. The dynamic power map is applied to the top layer (named the device layer hereafter) of the chip with the device layer thickness of 55.8 µm. For data collection of each building block, the thermal simulation is performed over the simulation domain that consists of the green blocks and the building block shown in Fig. 1. The simulation domain for data collection of each embedded building block is shown in Fig. 2. In such a setting, solution data collected from each building block is able to account for the variation of the block BCs induced by the power excitations outside the block.

TABLE I. TEMPERATURE INDEPENDENT MATERIAL PROPERTY. Specific heat, C Density, 𝝆 Thermal conductivity, k 751.1 (J/(kg•K)) 2330 (kg/m 3 ) 100 (W/(m•K))

All outer surfaces of the simulation domain for data collection are assumed adiabatic except for the bottom where the convection BC is implemented with a constant heat transfer coefficient and an ambient temperature 𝑇 456 of 45℃. The outer boundary surface of the simulation domain is 3 mm from the building block. The dynamic power density in each building block is randomly generated in time; in each time step, it represents an averaged power density over 48k CPU cycles at 3.5 GHz with a total power approximately equal to 8.9 W for L2 cache and 16 W for Core 1. In the surrounding of each building block, the dynamic power density is generated with a different random sequence, which offers the variation of the interface flux on each side of the building and allows the POD modes to adapt the realistic BC variation to construct a more effective POD model.

In this work, collection of thermal data in the CPU domain for training POD modes and calculations of POD model coefficients in ( 10)-( 13) are carried out in FEniCS-FEM. The solution of the ODE matrix equation in ( 9) and post-processing calculations for the predicted temperature in (1) are performed in C++ using solvers in the PETSc library.

A. Temperature data collection

In this work, L2 Cache and Core 1 of the quad-core chip are selected as Block 1 and Block 2, respectively, as shown in Fig. 1 and 2. The dynamic thermal simulation is performed for each of Block 1 and Block 2 in FEniCS-FEM independently to collect dynamic temperature data of each building block to generate its eigenvalues and POD modes. The eigenvalue represents the mean squared temperature variation captured by the corresponding POD mode, and therefore its spectrum reveals the information on the number of POD modes needed to offer accurate temperature solution. The eigenvalue spectrums of two building blocks are shown in Fig. 3. For both Block 1 and Block 2, a reduction in the eigenvalue by two orders of magnitude is observed from the first to the second mode and a decrease by four orders from the first to the third mode. Based on the rapid reduction in the eigenvalue for the first few modes, it is expected that the two-block POD model with a small number of modes is able to offer an accurate prediction of dynamic temperature solution. However, the expectation can be achieved only if the quality of the data collected from the DNSs is reasonably good.

B. Verification of the multi-block model

To demonstrate the validity of the two-block POD model, thermal simulation based on (9) for the domain consisting of Core 1 (Block 2) and its adjacent L2 Cache (Block 1), shown in Fig. 1, is performed. The dynamic power density applied to each of the 2 blocks is generated using a random sequence different from those used in the POD mode training. The adiabatic BCs are applied to all surfaces except for the bottom of the chip where a constant heat transfer coefficient is implemented with an ambient of 45 o C. Thermal simulation is also performed via FEniCS-FEM with identical settings including heat sources and BCs to validate the accuracy of the two-block POD model. In this demonstration, a same number of POD modes is used in both blocks.

The LS error estimated from the equation below for the twoblock POD model is a function of the number of POD modes.

where the index i denotes the time step (snapshot), and 𝑇 ! (𝑟 ⃑) and 𝑒 ! (𝑟 ⃑) are the temperature solution from FEniCS-FEM and the temperature difference between FEniCS-FEM and the POD model, respectively. For the two-block POD model, the DG method with an adjustable penalty number 𝑁 0 [28], [29] is used to enforce the thermal continuity across the interface between Core 1 and L2 Cache. The effect of penalty number on the LS error is shown in Fig. 4 for the two-block POD model. It is observed that, when using 2 or more modes in the POD model, the LS error reaches a minimum value with Nµ near 7. The LS error vs. the number of modes with Nµ = 7 is thus plotted in Fig. The POD simulation demonstrated above reveals a 4-order reduction in the numerical DoF, compared to FEniCS-FEM, which results in a significant saving in computing time. The POD simulation includes solving the ODE in ( 6) and the post processing calculation using (1) to recover the temperature solution. The computational time of thermal simulation for the selected two blocks using FEniCS-FEM and the two-block POD model is shown in Table Ⅱ, where Post1 and Post2 denote the post-processing calculations of temperature in the entire domain and device layer, respectively. As shown in Table Ⅱ, thermal simulation based on the two-block POD model with 3 modes is 1959 times faster than FEniCS-FEM. Practically, only the temperature in the device layer is required, which would offer a speedup of 3918 times, compared to FEniCS-FEM.

Based on the results presented in Fig. 4, the optimal penalty number is Nµ = 7, and thus the detailed comparison of the dynamic thermal distributions obtained from FEniCS-FEM and the two-block POD model is given below with Nµ = 7. As expected according to the eigenvalue spectrum shown in Fig. 3 and the LS error in Fig. 5, the temperature evolution predicted by the two-block POD model with just 3 POD modes is in very good agreement with that obtained from FEniCS-FEM over the entire simulation time. The temperature evolution in time at the center of L2 Cache is given in Fig. 6. Similarly, the dynamic temperature at the center of Core 1 is illustrated in Fig. 7, where the temperature solutions obtained from the two-block POD model with 3 or 5 POD modes and FEniCS-FEM almost overlap each other.

TABLE II. CONSUMPTIONAL TIME OF THERMAL SIMULATION FOR THE MULTI-BLOCK POD AND FENICS-FEM METHODS. The temperature distribution at t = 6.6 ms from L2 Cache to Core 1 along the centers of these 2 blocks across the interface is illustrated in Fig. 8. The temperature profile provided by the two-block POD model with 3 POD modes agrees very well with the temperature profile from FEniCS-FEM, which is consistent with the information indicated by the eigenvalue spectrum in Fig. 3 and the LS error in Fig. 5. Compared with the FEniCS-FEM results in the centers of L2 Cache and Core 1, approximately 2.6% and 1.0% (or 0.09 °C and 0.14 °C) differences, respectively, are achieved when using the two-block POD model with 3 POD modes. IV. CONCLUSION A multi-block thermal simulation methodology enabled by the data-driven POD model and domain decomposition has been investigated. The approach has been applied to develop a twoblock POD model for thermal simulation of 2 selected blocks, including Core 1 and L2 Cache, from a quad-core chip, AMD ATHLON II X4 610e. This study has shown that an appropriate penalty number Nμ is needed in the 2-block POD model to minimize the interface discontinuity for an optimal prediction of the dynamic thermal distribution in the 2-block domain. It is found that an LS error below 3% can be achieved for the 2-block POD model with 3 modes in each block if 4 ≤ 𝑁 0 ≤ 10. This results in a reduction in the numerical DoF by nearly 4 orders of magnitude and leads to nearly 2000 times or 4000 times of speedup for a thermal prediction of the entire 2-block domain or the device layer, respectively, compared to FEniCS-FEM. This work initiates the development of a multi-block POD thermal simulation methodology at the chip level for an entire CPU or GPU. Such a multi-block concept offers a very efficient approach to generation of the POD modes and calculations the POD model parameters to construct a POD model for very large chips like multicore CPUs or many-core GPUs, especially when a higher resolution is needed. All CPUs and GPUs are designed and constructed based on building blocks, such as cores, caches, I/O units, memory modules, etc., in the selected AMD ATHLON processor. When developing a multi-block POD model, a useful practice would be using these standard building blocks to partition the entire processor into a multi-block domain. For a semiconductor chip consisting of a large number of POD blocks, parallel computing can also be applied to further improve the POD simulation efficiency.