二维稳态热传导的问题CPU/GPU并行求解
【IT168 技术】
注：本文为IT168&NVIDIA联合举办的“如何并行化我的应用”方案征集活动第二名。本次方案征集活动详情见：http://cuda.itpub.net/thread-1299715-1-1.html。近期活动的大部分方案，将会逐步与大家分享，不可错过哦！
CUDA ZONE专区：http://cuda.it168.com/
CUDA技术论坛：http://cuda.itpub.net
热传导问题是材料成形模拟中经常需要求解的问题。热传导问题求解得到的是温度场。几乎所有材料成形模拟的过程都与温度场有关。对热传导问题的求解具有重要意义。
1.应用背景
制造业是一个国家的基础行业，它可以代表一个国家的工业水平。材料加工行业是制造业的重要组成部分。材料加工行业涵盖广泛，包括注塑、铸造、锻造、焊接、板料冲压等。现今，材料成形计算机模拟已经应用到材料加工行业的方方面面。材料成形模拟软件具有较强的实践指导意义，而且在规模不大问题上这种指导意义就可以体现出来。所以，国内的材料成形模拟大多仍停留在单机串行求解问题的水平上。但是，随着模拟工程的规模不断增长，工程应用对材料成形模拟软件的性能要求不断提高。为了满足实际需求，材料成形模拟软件需要实现并行化，以提高性能和求解大型问题的能力。
无论问题规模大小，几乎所有材料成形模拟的过程都与热传导问题有关。前面提到的注塑、铸造、锻造、焊接、板料冲压的模拟都与热传导问题密切相关。这类问题的求解过程中，热传导问题的求解总是不可忽略的步骤。实践发现，在对材料成形模拟问题进行求解时，目前的计算速度和客户理想需求的差距在20倍以上。
2.业务规模和瓶颈 以注塑模拟为例。
当网格规模为千万网格时，求解热传导问题需要80~200分钟，如果网格规模达到一亿，则需要5到10个小时。而客户的理想需求则分别为15分钟以内和60分钟以内。在存储需求上，千万网格的数据量大约是10GB，而一亿网格的数据量将超过100GB。对于千万网格的问题规模，网格数据可以存储在大型服务器的内存中，但对于一亿网格的数据量则必须存储在硬盘中。
根据以上描述的问题规模，热传导问题的计算瓶颈在于计算周期和数据存储。而并行计算，有望解决计算周期较长的问题。如果能对热传导问题的求解过程进行有效的并行化，则可以显著缩短材料成形模拟的整体周期。
3.并行策略
热传导问题的分类有多种。按是否与时间有关，可分为稳态问题和瞬态问题;按几何模型的维度，可分为二维和三维;按求解方法，又可分为Jacobi迭代法，fourier迭代法等。本文以具有代表性的二维稳态Jacobi迭代法求解热传导问题为例。
Jacobi方法是一种有限差分迭代方法。其迭代公式为
T(i, j) = T(i - 1, j) + T(i + 1, j) + T(i, j - 1) + T(i, j + 1)
在Jacobi方法的迭代中，迭代计算时仅仅用到上次迭代结果中的数据，当前各个网格处的迭代计算之间没有直接关系。所以Jacobi方法具有很好的并行性。串行计算时，所有网格处的迭代计算是串行执行的。Jacobi迭代的结束条件是，前后两次迭代的误差小于指定数值。
串行代码如下：
{
step++;
count = 0;
for(int i = 1;i < n - 1; i++)
for(int j = 1; j < n - 1; j++)
{
w[i][j] = (m[i - 1][j] + m[i + 1][j]
+ m[i][j - 1] + m[i][j + 1]) / 4.0;
if(fabs(w[i * n + j] - m[i * n + j]) < epsilon) count++;
}
temp = m; m = w; w = temp;
}
并行化的策略是尽可能并行化每一迭代步内所有网格处的计算。但有由于受到处理器个数和同步开销的限制，这样做事不现实的。较为实际的策略是，将所有网格分为合适数目的网格块，各个网格块的计算并行执行。
本文使用了MPI和CUDA两种方法分别进行并行计算。两种并行方法的并行策略稍有不同。MPI是适用于计算机集群的并行计算方法。在使用MPI实现并行计算时的策略如上所述。为了提高PI的并行效率，降低数据传输引起的时间损耗，具体实现时适用了异步传输函数。
{
step++;
//Synchronize boundary data
if (nodeIndex % 2 == 0)
{
//Next
if (nodeIndex != nodeCount - 1)
{
MPI_Isend(m + jobEnd * n - n, n, MPI_DOUBLE, nodeIndex + 1,
0, MPI_COMM_WORLD, &requestSendNext);
MPI_Irecv(m + jobEnd * n, n, MPI_DOUBLE, nodeIndex + 1,
1, MPI_COMM_WORLD, &requestReceiveNext);
}
//Previous
if (nodeIndex != 0)
{
MPI_Isend(m + jobStartingPoint * n, n, MPI_DOUBLE,
nodeIndex - 1, 1, MPI_COMM_WORLD, &requestSendPrevious)
MPI_Irecv(m + jobStartingPoint * n - n, n , MPI_DOUBLE,
nodeIndex - 1, 0, MPI_COMM_WORLD,
&requestReceivePrevious);
}
}
else
{
//Previous
MPI_Irecv(m + jobStartingPoint * n - n, n, MPI_DOUBLE,
nodeIndex - 1, 0, MPI_COMM_WORLD,
&requestReceivePrevious);
MPI_Isend(m + jobStartingPoint * n, n, MPI_DOUBLE,
nodeIndex - 1, 1, MPI_COMM_WORLD, &requestSendPrevious);
//Next
if (nodeIndex != nodeCount - 1)
{
MPI_Irecv(m + jobEnd * n, n, MPI_DOUBLE,
nodeIndex + 1, 1, MPI_COMM_WORLD, &requestSendNext);
MPI_Isend(m + jobEnd * n - n, n, MPI_DOUBLE,
nodeIndex + 1, 0, MPI_COMM_WORLD, &requestReceiveNext);
}
}
//Compute Inner Data
localCount = 0;
for(int i = jobStartingPoint + 1; i < jobEnd - 1; i++)
rowDataIteration(n, *step, epsilon, i, m, w, &localCount);
//Compute boundary data
int rowIndex;
if (nodeIndex == 0)
{
rowIndex = jobStartingPoint;
rowDataIteration(n, *step, epsilon, rowIndex, m, w, &localCount);
MPI_Wait(&requestReceiveNext, &statusNext);
rowIndex = jobEnd - 1;
rowDataIteration(n, *step, epsilon, rowIndex, m, w, &localCount);
MPI_Wait(&requestSendNext, &statusNext);
}
else if (nodeIndex == nodeCount - 1)
{
rowIndex = jobEnd - 1;
rowDataIteration(n, *step, epsilon, rowIndex, m, w, &localCount);
MPI_Wait(&requestReceivePrevious, &statusPrevious);
rowIndex = jobStartingPoint;
rowDataIteration(n, *step, epsilon, rowIndex, m, w, &localCount);
MPI_Wait(&requestSendPrevious, &statusPrevious);
}
else
{
//还可以MPI_Waitany()来提高性能***
MPI_Wait(&requestReceivePrevious, &statusPrevious);
rowIndex = jobStartingPoint;
rowDataIteration(n, *step, epsilon, rowIndex, m, w, &localCount);
MPI_Wait(&requestReceiveNext, &statusNext);
rowIndex = jobEnd - 1;
rowDataIteration(n, *step, epsilon, rowIndex, m, w, &localCount);
MPI_Wait(&requestSendPrevious, &statusPrevious);
MPI_Wait(&requestSendNext, &statusNext);
}
totalCount = 0;
MPI_Allreduce(&localCount, &totalCount, 1,
MPI_INT, MPI_SUM, MPI_COMM_WORLD);
temp = m; m = w; w = temp;
}
CUDA是NVIDIA公司推出的基于GPU的通用计算工具。GPU并行计算是利用GPU中的简化计算核心SP(流处理器)来进行计算，由于单个GPU中集成了大量SP，CUDA通过利用这些SP协同计算，实现较高的加速比。由于CUDA会自行将计算线程映射到SP，而且线程越多CUDA的并行效率越高。所以基于GPU的Jacobi并行采取了完全并行，即同一迭代步内，对所有网格的计算都事实并行。另外，对于迭代的结束条件也进行了并行优化。在串行算法中，迭代结束条件即是对前后两次迭代结果的差值求最大值。在GPU并行时，采用了并行的求最值算法，使得收敛判断(即迭代结束条件)也得到了很好的并行。
GPU(CUDA)代码如下：
__global__ void kernelJacobiIteration(const int n, double *m, double *w)
{
int id = blockIdx.x * blockDim.x +threadIdx.x;
if (id < (n - 2) * (n - 2))
{
int row = id / (n - 2);
int column = id - row * (n - 2);
int location = (row + 1) * n + (column + 1);
w[location] = (m[location - 1] + m[location - n]
+ m[location + 1] + m[location + n]) / 4.0;
}
}
//收敛判断CUDA代码
//kernel of getting epsilon between d_m and d_w
__global__ void kernelGetEpsilon(const int n,
double *m, double *w,
double *ep)
{
int id = blockIdx.x * blockDim.x +threadIdx.x;
if (id < (n - 2) * (n - 2))
{
int row = id / (n - 2);
int column = id - row * (n - 2);
int location = (row + 1) * n + (column + 1);
ep[id] = fabs(m[location] - w[location]);
}
}
//kernel of max of a num group
__global__ void kernelGetMax(const int count, double *ep)
{
int id = blockIdx.x * blockDim.x + threadIdx.x;
if (id < count)
if (ep[id] < ep[id + count]) ep[id] = ep[id + count];
}
4.并行效果
并行测试数据如下：
测试数据中，在网格规模为1.6千万时，串行计算耗时2741秒(约45分钟)，MPI并行(2个节点)耗时1427秒(约24分钟)，GPU并行(NVIDIA GTX260+，216个SP)耗时171秒(约3分钟);MPI和GPU并行加速比分别是1.92和16.03。
总的来说，热传导问题的并行求解取得了很好的效果。虽然该算法较为简单，但表明在热传导问题上，并行求解将会有较好的效果。
