Solving the inverse heat conduction problem using NVLink capable Power architecture

Direct heat conduction problem

There are various fundamental modes of heat transfer, but this paper deals only with transient conduction when the temperature of the workpiece changes as a function of time. The determination of the temperature of any points of the object at any moment often calls for some computer-assisted approximation methods. In the case of three-dimensional objects, these calculations can be very complicated and resource intensive, preventing the efficient usage as a GA cost function. For cylindrical objects, an essential simplification can significantly decrease the computational effort. As can be seen in Fig. 1 it is enough to model the middle cross-section of the cylinder, resulting in the two-dimensional axis-symmetrical heat conduction model. The radius of the cylinder is noted by R and Z.

The cylinder is subjected to a longitudinal local coordinate and time varying Heat Transfer Coefficient HTC(z, t) on all its surfaces. Both the thermal conductivity, density and the heat capacity are varying with the temperature, k(T), ρ(T) and C_p(T). It has to be noted that phase transformations of the materials applied do not occur during the experiments, therefore latent heat generation induced by phase transformations is not considered.

Based on this simplified model, the mathematical formulation of the nonlinear transient heat conduction can be described as Eq. (1), with the following initial and boundaryconditions Eqs. (2)–(5):

(1)${\displaystyle \frac{\partial }{\partial r}\left(k\frac{\partial T}{\partial r}\right)+\frac{k}{r}\frac{\partial T}{\partial r}+\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)+q_v=\rho C_p\frac{\partial T}{\partial t}}$ (2)${\displaystyle T\left(r=0,z,t=0\right)=T_0,}$ (3)${\displaystyle k\frac{\partial T}{\partial z}{|}_{0\le z\le L,r=R}=HTC\left(z,t\right)\left[T_q-T\left(r=R,z,t\right)\right],}$ (4)${\displaystyle k\frac{\partial T}{\partial z}{|}_{z=0,0\le r<R}=HTC\left(z=0,t\right)\left[T_q-T\left(r,z=0,t\right)\right],}$ (5)${\displaystyle k\frac{\partial T}{\partial z}{|}_{z=Z,0\le r<R}=HTC\left(z=Z,t\right)\left[T_q-T\left(r,z=Z,t\right)\right],}$ where

• r, z—local coordinates;

• t—time;

• R—radius of the workpiece;

• ρ(T)—density of the object;

• T₀—initial temperature of the workpiece;

• T_q—temperature of the cooling medium;

• T(r, z, t)—temperature of the workpiece at given location/time;

• k(T)—thermal conductivity (varying with temperature);

• C_p(T)—heat capacity (varying with temperature);

• HTC(z, t)—heat conduction (varying with local coordinate and time).

Discretising of these equations using the weighted Schmidt explicit finite difference method, results in nine different equations, according to the location within the object.

• case (a)—inner items surrounded by four neighbouring items;

• case (b)—inner items at the centre line of the object;

• case (c)—boundary items at the outer surface of the object;

• case (d)—boundary items at the top surface of the object;

• case (e)—boundary items at the bottom surface of the object;

• case (f)—boundary item at the outer top corner;

• case (g)—boundary item at the inner top corner;

• case (h)—boundary item at the outer bottom corner;

• case (i)—boundary item at the inner bottom corner.

The reliability of the estimated HTC strongly depends on the numerical approach applied. During the DHCP calculations, the resolution of the finite grid is essential, parameters n = 10, m = 34 were used, where n is the number of points horizontally, and m is the number of points vertically. A sufficiently small dt (simulation time interval) value is also necessary to ensure the accuracy of the method (dt = 0.01 second).

Massively parallel solution

It is necessary to solve the heat transfer equations for each finite item for the consequent time periods. Equations for calculating the heat movement between neighbouring items within the same time can be solved in a parallel fashion because none of these calculations modifies the output of the others. Moreover, for several items, the steps of these calculations are the same, only the input parameters are different. This makes it possible to implement an efficient data-parallel algorithm. Assuming n × m number of threads, the parallelization level is:

• (n − 2) × (m − 2) threads have to calculate case a;

• (m − 2) threads have to calculate case b;

• (m − 2) threads have to calculate case c;

• (n − 2) threads have to calculate case d;

• (n − 2) threads have to calculate case e;

• 1 + 1 + 1 + 1 threads have to calculate case f, g, h and i.

Based on this data parallel fashion, it is worth implementing a GPU-based implementation. Every GPU thread is responsible for one item in the finite grid. The task of these threads is the determination of the temperature in the corresponding location in each consecutive time step. These threads can work in parallel with the equations of the same time step. Nevertheless, these are not totally independent, because as an input parameter, threads need the temperature of the neighbouring elements at the previous time step. The threads have to wait until all the others have finished the calculations for the actual period to fulfil this before they can continue working on the next time step.

This needs an explicit synchronisation after every time step, and in the CUDA framework, the only efficient barrier synchronisation is given by the __syncthreads() function, which synchronises the threads within the same block. According to this, it is necessary to keep all threads working on the same finite grid in the same CUDA block. There is an upper limit for the number of threads within a block (1,024 for current architectures); but, in the case of ordinary configuration, this is not a limiting factor (10⋅34 = 340 threads were used).

Running 340 threads is not sufficient to fully utilise a modern GPU. One P100 accelerator card has 3,584 CUDA cores. Therefore, this low number of threads leads only a low theoretical occupancy level (number of used cores/number of available cores = 340/3,584 = 9.48%), and the practical utilisation is expected to be even worse (not to mention, that one server node has four cards installed). At this point, recall that this DHCP calculation is responsible for the cost calculation part of a genetic algorithm. A population containing P chromosomes needs P number of thermal history generations in the evaluation phase of the GA. These are independent calculations, so it is feasible to launch these in a parallel fashion using more than one multi-processors. Given this higher level of parallelism, the required number of threads becomes P × 340, which is enough to utilise the computing power of the graphics accelerators fully.

This observation is the key point for the multi-GPU implementation. The already mentioned independence makes it possible to distribute these fitness calculations among the available devices. In the case of G number of GPUs, it is worth to assign each chromosome to one of the GPUs using Eq. (6), which shows the number of assigned chromosomes (P_i) for the ith GPU (where 1 ≤ i ≤ G). (6)${\displaystyle P_i=\left\{\begin{array}{cc}P-\left(G-1\right)⌊\frac{P}{G}⌋,\hfill & \text{if}i=1\hfill \\ ⌊\frac{P}{G}⌋,\hfill & \text{otherwise}\hfill \end{array}\right.}$

Further optimisation

Warp divergence

When a thread block is assigned to a multiprocessor, it is divided into further groups of 32 threads called warps. This partitioning is predetermined: threads with index 0–31 will be the first warp, threads in the 32–63 index range the second warp, and so on NVIDIA (2014). The warp is the lowest unit of execution scheduling, at any time, one instruction is fetched and executed by all threads within the same warp. The operands of these instructions can be different; therefore, in the case of conditional statements, different threads in the same warp should take different paths. The solution for this warp divergence is that all threads within a warp must take both branches of the conditional statement; threads with false condition are stalling, while the others are executing the instructions of the true branch, and threads with true condition are stalling, while the others are executing the false branch. These stall phases can significantly decrease the occupancy; hence, warp divergence should be avoided to obtain the best performance.

In the case of DHCP calculations, the number of threads is equal to the size of the finite grid (n × m). Every thread corresponds to the calculations of one item in the grid, and the first intuition shows that thread indices should be the same as the indices of the corresponding finite item. However, this leads to heavy warp divergence; for example, threads of the first warp (identified by indices below) should execute the following code paths (Fig. 2):

• (0, 0) thread have to calculate case g;

• (9, 0) thread have to calculate case f;

• (1, 0)–(8, 0) threads have to calculate case d;

• (0, 1)–(0, 3) threads have to calculate case b;

• (9, 1)–(9, 2) threads have to calculate case c;

• (1, 1)–(8, 2) and (1, 3) threads have to calculate case a.

Separation of thread indices and the corresponding finite grid locations makes it possible to decrease this divergence significantly. Threads inside the first warp should calculate the equations corresponding to the elements at the centre line (32 items). The second warp should correspond to the outer surface elements (32 items). Threads in the next eight warps should be responsible for the inner elements (272 elements), and the last warp should solve all the remaining equations (20 elements).

In this partitioning, only the last warp has divergent threads. The tests show that the gained speed-up is about 10–20%.

 
________________________________________________________________________________________________________ 
Algorithm 1 Data parallel DHCP solver 
________________________________________________________________________________________________________ 
Require: HTC: Heat Transfer Coefficient 
Require: T0: initial temperature 
Require: R: reference point inside the work object 
Ensure: S[ ]: recorded temperature values at the reference point 
 1: function DHCP-kernel( )                  ⊳ threadId: unique identifier in 0...n ⋅ m − 1 
 2:   (i,j) ← AssignThreadToFiniteNode(threadId)                  ⊳ ‘Warp divergence’ 
 3:   Ti,j ← T0         ⊳ Shared array to store actual temp values (‘Using shared memory’) 
 4:   synchronise threads 
 5:   for t ← 0 to S.length − 1 do 
 6:        time ← t ∗ dt                                                          ⊳ Actual time 
 7:        temp ← Ti,j                         ⊳ Thread-level variable (actual temp at (i,j) pos) 
 8:        switch (i,j) do ⊳ Calculate heat movement (‘Direct heat conduction problem’) 
 9:              case (0,0) 
10:                   temp ← temp + HeatTransferInnerTopCorner(T,HTC) 
11:              case ... 
12:                   ... 
13:        end switch 
14:        synchronise threads 
15:        Ti,j ← temp 
16:        synchronise threads 
17:        if threadId = 0 then                         ⊳ First thread calculates the result 
18:              S[t] ← InterpolateTempAt(T,R)                 ⊳ Calculate temp at ref. point 
19:        end if 
20:   end for 
21: end function 
22: 
23: function DHCP(HTC,T0,R) 
24:      AllocateGPUMemory( ) 
25:      CopyFromHostToDevice(HTC,T0,R) 
26:      DHCP-Kernel ⋘ n ∗ m ⋙ ( )                        ⊳ Launch n ⋅ m parallel threads 
27:   CopyFromDeviceToHost(S) 
28:   FreeGPUMemory( ) 
29:   return S 
30: end function 
________________________________________________________________________________________________________    

Benchmarking methodology

Several benchmarks were run with the CPU and the GPU implementations focusing on the following questions:

• What is the correlation between the number of CPU cores/GPU devices and the required runtime? A linear correlation is expected because of the weak dependencies between the different tasks.

• In the case of GAs, which hardware configuration is preferred for a given population size? The expectation is that it is worth using the CPU for small populations and the GPU for larger populations.

• The amount of input parameters (HTC control points) and output results (fitness values) is relatively small compared to other HPC applications. Is the new NVLink technology between the CPU and GPU able to significantly reduce the memory transmission time for these small data transfers?

The details of the test environments are as follows:

• Test Environment 1 (TE1)

  • CPU—1 × Intel(R) Core(TM) i7-2600

    • 4 physical cores

    • 8 logical cores

  • GPU—2 × GeForce GTX TITAN Black

    • CUDA cores: 2880

    • Memory: 6 GB

    • Link: PCIe 4× + PCIe 8×

• Test Environment 2 (TE1)

  • CPU—2 × IBM POWER8NVL

    • 10 cores

  • GPU—4 × Tesla P100-SXM2

    • CUDA cores: 3584

    • Memory: 16 GB

    • Link: PCIe 8×+ NVLink

TE1 was a Windows 10 based desktop machine with two GeForce GTX TITAN Black cards installed into the PCIe slots. The IHCP codebase was developed using Visual Studio 2015 (with CUDA 8.0) and compiled with the Nvidia nvcc compiler to standard 64 bit executables. These binaries were launched from the standard Windows command prompt.

The C + + standard std::chrono:high_resolution_clock object was used to measure the execution time. To decrease uncertainty, 20 independent tests were run for all parameter sets, removing the lowest and highest runtimes (5–5%) and computing the average of the remaining values.

TE2 was an Ubuntu Linux based Power System S822LC node in the IBM’s Accelerator Lab for High Performance Computing cluster. This is a two socket system equipped with the followings: 2 POWER8 10 Core CPUs, up to 1TB System Memory, and 4 NVIDIA P100 GPUs with NVLink connection. This system includes the exclusive NVLink high bandwidth interconnect between the POWER8 CPUs and NVIDIA GPUs, providing unprecedented bandwith between the CPU and GPUs not available on any other platform. Each P100 GPU is connected to one of the CPUs with 80 GB/sec bi-directional bandwith NVLink connections, and each pair of P100 GPUs are also directly connected to each other with 80GB/sec bi-directional bandwidth NVLink connections (Fig. 3). The same codebase was compiled with gcc and CUDA 8.0 tools using the arch=compute_60,code=sm_60 flags. The benchmarking methodology was the same as for TE1.

Several tests were run using different population sizes (P), where P = 10, 20, 30…2,000. As a second parameter, the number of GPU devices (G) was changed, too. G = 1, 2 for the first environment, and G = 1, 2, 3, 4 for the second environment. To compare the CPU and GPU performance, the same tests were run using a different number of CPUs (C), where C = 5, 10, 20. The executed steps of the DHCP process are not affected by the actual HTC values. Therefore, predetermined input values were used for the simulations (instead of random chromosomes from a real genetic algorithm).

GPU runtimes

Table 1 and Fig. 4 show the GPU results for TE1, and Table 2 and Fig. 5 show the GPU results for TE2.

As expected, the required runtime inversely linearly correlates with the number of devices. However, in practice, this is a bit more complex, because the series showing the runtimes in both figures are not straight lines, which is caused by the speciality of GPU hardware. This fragmentation is best visible in the case of the one GTX Titan Black card (black solid series in Fig. 4). The total runtime is slightly increased from population sizes 10 to 60, but the required runtime for population size 70 is almost twice as large. The explanation for this phenomenon is based on the GPU architecture. The GTX Titan Black has 15 streaming multiprocessors (each of them containing 192 processing units), and every multiprocessor can execute four heat transfer simulations simultaneously. Accordingly, launching one or 60 parallel simulations takes a similar time. In the case of a greater number of threads, the scheduling (including memory transfers) becomes more complex; consequently, the effect of these steps becomes less sharp.

As visible in Fig. 5, these steps already exist for one P100 card, the runtime increases in every 50–60th step. The P100 has 56 multiprocessors, each of these can execute one simulation at the same time; therefore, the device can run 56 thermal history generations in parallel. This regularity can also be observed in the case of multiple GPUs; for example, in the case of 4 P100 devices, steps from 220 to 230 (56⋅4 = 226) and from 440 to 450 (56⋅4⋅2 = 448) has a significant impact to the runtime. It is worth evaluating as large a population as possible at the same time. Thus, the recommended population size should be near to these limits from below (220 or 440).

CPU runtimes

Table 3 and Fig. 6 show the CPU results for TE2.

CPU performance analysis is not in the focus of this paper, but these benchmarks have been run only for the CPU–GPU comparison. As visible, the increase in the number of cores effectively increases the performance (each core is responsible for one heat transfer simulation). On the other side, as the population size was increased, the runtime increased almost linearly.

In the case of CPU implementations, the implementation is much simpler. There is no need to transfer input data from the host to the device and the output data from the device to the host, and kernel launch overhead is also missing. Therefore, the expectation was that the CPU would be faster in the case of small population sizes (where the GPUs cannot take advantage of the high number of multiprocessors), but as is visible, this was not true.

Comparing these results to the GPU runtimes, it is evident that it is not worth using the CPU for any population size. In the case of 20 parallel heat transfer simulations (which is the ideal configuration for the server with 20 CPU cores), all P100 GPU implementations were 8–9 times faster. Above this population size, the difference becomes even bigger. In the case of 2,000 parallel simulations, the P100 cards were 123 times faster.

It also raises the question of whether it is worth implementing and using a hybrid solution combining the CPUs and GPUs together. For small population sizes, the answer is no. For example, in the case of 200 chromosomes, if the GPUs evaluate 199 of them (needs about 134 ms) and the CPUs evaluate only one (needs about 1,024 ms), the overall runtime will be higher (max(134 ms, 1,024 ms) = 1,024 ms) than the plain GPU implementation. The only occasion when it is worth considering the hybrid implementation is a large population size near the previously explained runtime steps from above. For example, in the case of 1,690 chromosomes and two GPUs, it would be worth to assign 20 of them to the CPUs. The CPU runtime will be about 1,024 ms, the GPU runtime for the remaining is 1,054 ms. Consequently, the overall execution time is 1,054 ms. Using the GPU implementation exclusively, it requires more, 1,132 ms. In the case of four GPUs, it would also be possible to find a similar situation, but not in the examined 10–2,000 population size interval.

Data transfer rates

During the tests, the data transfer rates were also recorded in the case of TE1 and TE2. Tables 4–5 and Figs. 7–8 show these results. The advance of the new NVLink based Power CPU–GPU connection is the high transfer bandwidth between the CPU and the GPU memory. As visible from the results, this works well in practice.

Table 4 shows that the measured results for TE1 were far from the theoretical maximal transfer speed and it is also visible, that using two GPUs have advantages compared to the single GPU configuration. These devices are installed into different PCI-E slots, so their data transfers can be processed in parallel (both cards were installed into PCI-E 16 × slots, but one of them works at 4 × speed, and the other one at 8 × speed). In the case of the CPU–GPU NVLink architecture (Table 5), the bandwidth was significantly higher (however, the amount of data is still low to achieve the theoretical maximum).

It is worth noting that with the two GeForce GPUs, the memory transfer time is well over twice as fast as with one card; in sharp contrast, for two P100 cards it is only about 40% faster than that with one card. And what is more interesting, the usage of three or four devices does not have any benificial effect. The reason for this is that the used P100 node has a special architecture in that one CPU socket and a maximum of two GPUs form an “island” in which any pair of the triad has fast data transfer. In contrast, transfer of data between islands is slower. In the case of N_G number of GPUs, the IHCP solver code uses N_G individual CPU threads to manage the data transfers between the host and the devices. This practice was satisfactory for TE1, but leads to the experienced strange behaviour for TE2, because all of these threads were scheduled on the same physical processor. To reach the maximum transfer rate with four devices, multiple CPU threads should be launched and bounded to logical CPUs that are appropriately placed for the GPUs we intend to use. However, customizing the already existing codebase to a specified architecture is out of the scope of this paper, especially in view of the fact that the NVLink capable server was even significantly faster without any fine-tuning.

Comparing the GTX and P100 based configurations, the latter was 12–49 times faster. The new NVLink architecture is, therefore, faster, and its overall transfer rate becomes higher with multiple cards. This speed up is achievable in the case of small population sizes (it is 16 times quicker in the case of 20 chromosomes and 2 GPUs). Therefore, this novel architecture is highly recommended for running the IHCP solver GAs.