CN101727653A - 一种基于图形处理器的多组分系统离散模拟计算方法 - Google Patents

Abstract

本发明提供一种基于图形处理器的多组分系统离散模拟计算方法，该方法在由多个计算节点组成的并行计算系统上实现，计算节点上安装有图形处理器GPU；包括：对所要模拟的多组分系统做空间分解；建立多组分系统中的粒子与网格之间的映射关系；将粒子信息分配给多个计算节点上的GPU内的线程块，并将粒子的作用势参数、模拟设置信息保存到GPU的内存中；在所述GPU上计算粒子所受到的作用力，并根据所述作用力对粒子的速度与位置进行更新；在对所述多组分系统模拟结束前，重复在GPU内完成的执行步骤，直到模拟结束后将模拟结果从所述并行计算系统中输出。本发明克服了GPU在计算多种作用势时由于选择导致的低效问题，与CPU相比，在计算效率上有很大的提高。

Description

一种基于图形处理器的多组分系统离散模拟计算方法

技术领域

本发明涉及多组分系统离散模拟的计算，特别涉及到在图形处理器GPU上实现多组分系统离散模拟计算的相关方法。

背景技术

多组分系统是指由两种以上不同性质的组分所形成的系统。在对多组分系统做模拟时，多组分系统中的不同组分可以用不同种类的粒子来表示，所述的模拟过程需要处理粒子运动以及粒子间的相互作用。在模拟过程中，粒子间的相互作用一般可以用二体作用“势”或者作用力函数表示。如在分子动力学MD模拟中，常见的Lennard-Jones势采用下列公式表示：

$\mathrm{\φ}\left(r_{\mathrm{ij}}\right)=4\mathrm{\ϵ}[a{\left(\frac{\mathrm{\σ}}{r_{\mathrm{ij}}}\right)}^{12}-b{\left(\frac{\mathrm{\σ}}{r_{\mathrm{ij}}}\right)}^6]$

其中，ε为作用势能，σ为作用特征长度(一般为分子直径)，r_ij为分子i到分子j的距离；系数a、b表示斥力项和引力项的相对大小，一般情况a＝b＝1，而a＝1和b＝0时表示只有斥力作用。

而MD模拟中的作用力函数则可以用下列公式表示：

$F_{\mathrm{ij}}=-\▿\mathrm{\φ}\left(r_{\mathrm{ij}}\right)=\left(\frac{48\mathrm{\ϵ}}{{\mathrm{\σ}}^2}\right)[a{\left(\frac{\mathrm{\σ}}{r_{\mathrm{ij}}}\right)}^{14}-\frac{b}{2}{\left(\frac{\mathrm{\σ}}{r_{\mathrm{ij}}}\right)}^8](r_i-r_j)$

其中，F_ij表示分子i受到分子j的作用力。

关于分子动力学MD模拟的具体情况可参见参考文献1：“Frenkel Daanand Smit Berend，1996.Understanding Molecular Simulation：FromAlgorithms to Applications.Orlando，Academic Press”。

现有技术中，用计算机实现前述的对多组分系统的模拟，主要是在计算机上实现了对粒子间势函数或者作用力函数的计算。在实际模拟过程中，由于所要模拟的粒子数目十分庞大，因此对计算机的计算能力的要求很高。现有技术中通常采用由中央处理器(CPU)为计算部件的计算机实现对多组分系统的模拟。但由于CPU自身性能在发展上的瓶颈，使得此类计算机并不能够很好地满足对多组分系统模拟的需要，带来计算效率较低，所能模拟的多组分系统的规模偏小等一系列问题。与CPU相比，图形处理器GPU集成了更多晶体管，而且有更高比例的晶体管用于数值计算，因此其计算速度远高于主流CPU，两者在计算能力上的差别在10～100倍左右，目前最先进的图形处理器GPU的计算能力已经达到了1Teraflops(万亿次)。因此，GPU特别适合并行度高的计算密集型计算领域，即采用GPU的计算机较采用CPU的计算机更适合多组分系统的计算机模拟。

现有技术中存在多种在计算机上实现多组分系统模拟的相关算法，但这些算法都是在以CPU为计算部件的计算机上实现的，由于CPU与GPU在结构以及工作原理等方面的诸多差异，使得这些算法并不能够直接移植到以GPU为计算部件的计算机上，即便能够移植，也会对模拟效果产生影响。例如：

1、在现有算法中，不同作用势的选择主要通过各个条件选项分支完成，但该方式对于流水线的优化极为不利。由于图形处理器GPU擅长顺序执行，但缺乏逻辑判断的硬件实现，GPU需完整计算作用势选择的各个分支后根据条件判断选择合适的分支，如此对于存在较多作用势的多组分系统，GPU会将大部分计算消耗在无用分支上，GPU计算发挥的性能将大打折扣。

2、由于GPU在硬件上缺乏跳转(jump)指令实现，因此目前的GPU计算还不能实现由函数指针来选择不同的作用函数，这也为现有算法到以GPU为计算部件的平台上的移植带来了困难。

发明内容

本发明的目的是克服现有的多组分系统模拟计算方法无法直接移植到GPU平台上的缺陷，从而提供一种在带有GPU的并行计算系统上实现计算的相关方法。

为了实现上述目的，本发明提供了一种基于图形处理器的多组分系统离散模拟计算方法，该方法在由多个计算节点组成的并行计算系统上实现，所述计算节点上安装有图形处理器GPU；该方法包括：

步骤1)、对所要模拟的多组分系统做空间分解，在所述空间分解中，将所述多组分系统的空间分为多个子区域，并将所述子区域分为网格；

步骤2)、建立所述多组分系统中的粒子与所述步骤1)所得到的网格之间的映射关系；

步骤3)、将所述多组分系统中的粒子信息按照所述空间分解的结果分别分配给所述多个计算节点上的GPU内的线程块，并将粒子的作用势参数、模拟设置信息保存到所述GPU的内存中；

步骤4)、在所述GPU上计算粒子所受到的作用力，并根据所述作用力对粒子的速度与位置进行更新；

步骤5)、重复对粒子的上述计算过程，直到对所述多组分系统模拟结束后将模拟结果从所述并行计算系统中输出。

上述技术方案中，所述的步骤4)包括：

步骤4-1)、在所述GPU的各个线程块上并行计算各个粒子的当前速度与当前位置，根据计算结果实现对粒子位移以及粒子速度的更新；

步骤4-2)、根据步骤4-1)的更新结果在各个GPU内实现所述粒子到所述网格映射的更新；

步骤4-3)、利用所述计算节点上的CPU进程实现GPU间的信息传递，实现粒子在不同GPU间的迁移以及边界粒子传递；

步骤4-4)、根据所述粒子的类型，选择相应的粒子作用势参数，实现粒子间作用力的计算；

步骤4-5)、根据所述粒子间作用力的计算结果，再次对所述粒子的速度做更新操作。

上述技术方案中，在所述的步骤1)中，按照所述并行计算系统中的GPU个数实现对所述子区域的划分，使得所述子区域的个数与所述GPU的个数相等；在所述子区域内按照作用截断距离信息划分网格。

上述技术方案中，在所述的步骤3)中，将所述多组分系统中的粒子信息按照所述空间分解的结果分别分配给所述多个计算节点上的GPU内的线程块时，将所述多组分系统的一个子区域分配给一个所述的GPU，将所述子区域内的一个网格分配给所述GPU内的一个线程块，所述网格内的一个粒子由所述线程块内的一个线程做专门的处理。

上述技术方案中，在所述的步骤3)中，所述的作用势参数包括作用截断距离r、作用特征长度σ、作用势能ε、用于表示斥力项相对大小的系数a、用于表示引力项相对大小的系数b在内的多个参数。

上述技术方案中，所述的步骤4-2)包括：

步骤4-2-1)、根据粒子的新位置计算粒子在运动后所在网格的坐标；

步骤4-2-2)、各个网格备份原有的粒子到网格的映射信息，并从所述映射信息中找出网格坐标与当前网格坐标不相符合的粒子，删除该粒子；

步骤4-2-3)、网格根据周边邻居网格所备份的粒子到网格的映射信息，从周边邻居中找出属于自己的粒子，并将这些粒子放入本网格的粒子映射信息中。

上述技术方案中，所述的步骤4-3)包括：

步骤4-3-1)、GPU线程块中的各个线程分别将位于内边界区域和外边界区域中的各个粒子的到网格的映射信息保存到所述GPU内存的缓冲区内；

步骤4-3-2)、将GPU内存缓冲区内的信息复制到对所述GPU进行管理的CPU进程所对应的发送缓冲区内；

步骤4-3-3)、由所述CPU进程将所述发送缓冲区内的信息传输到目标邻居进程的接收缓冲区，将其他邻居进程传入的信息放置于本进程的接收缓冲区内；

步骤4-3-4)、将所述接收缓冲区内的信息复制到所述GPU的内存缓冲区；

步骤4-3-5)、GPU解析内存缓冲区内的粒子信息，将跨界粒子添加到内边界网格粒子的映射信息中，或者对位于外边界网格的粒子的映射信息做替换。

上述技术方案中，在所述的步骤4-3-3)中，CPU进程与邻居进程间的通信采用shift模式。

上述技术方案中，所述的步骤4-4)包括：

步骤4-4-1)、GPU在内核中为一个线程块开辟共享内存；

步骤4-4-2)、在所述的共享内存内存放该线程块所对应的网格的信息以及该网格的邻居网格的信息；

步骤4-4-3)、在共享内存内，根据粒子类型，选择相应的粒子作用势参数，计算线程块所对应网格内的粒子与一个邻居网格内粒子的所有作用力；

步骤4-4-4)、在计算完一个线程块所对应网格内的粒子与一个邻居网格内的粒子间的相互作用力后轮换邻居网格，重复上述的步骤4-4-3)，直到执行完与周边所有邻居粒子间的作用力计算后才将线程块对应网格内的包括作用力大小的粒子信息写回。

本发明的优点在于：

1、本发明的方法预先将作用势参数进行存储，然后在计算粒子对相互作用力时，根据涉及粒子的种类属性，通过简单数值计算得到所需作用势参数在存储时的位置，从而根据所获得的作用势参数完成粒子间的作用计算，克服了GPU在计算多种作用势时由于选择导致的低效问题。

2、本发明的方法在多个GPU的多个线程块内并行实现对所要模拟的多组分系统中的粒子的计算，与CPU相比，在计算效率上有很大的提高。

附图说明

图1为在本发明的一个实施例中所提供的并行计算系统的示意图；

图2为本发明的方法的流程图；

图3为在本发明的方法中实现线程块与网格对应的示意图；

图4为本发明的方法中对并行计算区域进行分解的示意图；

图5为本发明在一个实施例中在邻居进程间所采用的shift模式的示意图；

图6为采用本发明的方法对气液固三相体系进行模拟的结果示意图；其中图6(a)表示在t为0时的结果，图6(b)表示t在300时的结果，图6(c)表示t在500时的结果，图6(d)表示t在900时的结果，图6(e)表示t在1250时的结果。

具体实施方式

下面结合附图和具体实施方式，对本发明做进一步说明。

本发明在模拟多组分系统时，以GPU为计算部件的计算机系统为基础，实现模拟过程。虽然单个的GPU与对应的CPU相比，在计算性能上已经有了很大的提高，但为了更好地对多组分系统中大量粒子的模拟，本发明采用由多个GPU组成的并行计算系统实现对多组分系统的模拟计算。所述的并行计算系统由多个计算节点组成，在每个计算节点上都分别安装有所述的图形处理器GPU，不同计算节点间通过网络连接。

在本发明的一个实施例中，提供了如图1所示的并行计算系统，该计算系统包括4个计算节点，每个计算节点安装有2块

公司生产的TeslaTM C870图形处理器(GPU)，由Linux系统构建计算环境，采用以太网连接各计算节点，其中的GPU包括计算内核、多处理器片下内存以及多处理器片上内存。多处理器片下内存则包括全局内存、常量内存、纹理内存在内的多种内存，多处理器片上内存包括寄存器、当地内存以及共享内存。多处理器片下内存中的全局内存和常量内存也被称为GPU显存。

在本实施例中，假设用上述的并行计算系统对一种三组分系统做分子动力学MD模拟，所要模拟区域的大小为80×160×80，在该模拟区域内，具有厚度为3的两边壁面(由76800个数密度n_s＝1的固体分子组成)、固定在上部的直径D_s＝20的固体颗粒(含3544个固体分子)、位于底端的直径D_b＝50的气泡(含45965个气体分子)以及填充在其他区域内的数密度n₁＝0.733的液体分子(总数N₁＝643240)。上述三组分系统的初始温度设定为k_BT＝1。

下面参考图2，对如何在图1所示的并行计算系统上实现对三组分系统的模拟进行说明。

将所要模拟的三组分系统的上述信息读入到前述的并行计算系统中，在所述的计算机系统中实现对模拟区域的空间分解。从前面的说明中已经知道，本发明为了提高效率，将多个具有GPU的计算节点通过网络并行连接，而各个GPU中的线程间又具有并行计算的功能。为了提高模拟计算的效率，需要将模拟区域中的不同空间分配给不同的GPU以及GPU中的不同线程进行相关的计算，这种分配也被称为任务划分。因此，首先需要对模拟区域做空间分解，即将三维的模拟区域做立体分割，得到相应的子区域。在本实施例中，采用二维分割的方式将所述的模拟区域分割成Nx×Ny个子区域，在各个子区域内再根据作用截断距离信息(作用截断距离是指当两个粒子间距离大于该作用截断距离时，则两粒子之间就被视作不发生作用)划分网格。在划分子区域时，所划分子区域的个数通常与并行计算机系统内的GPU个数相关，一个GPU对应一个子区域。本实施例中，由子区域进一步划分所得到的网格的大小为5(此处的大小为无量纲值，可以根据实际体系需要转化为带单位的真实值)，根据粒子数密度(一般在0.6～1.1之间)等信息可知大小为5的网格内粒子数在100个左右。

在对模拟区域做空间分解后，并行计算机系统根据模拟区域中的粒子所在的位置，在粒子与粒子所在网格间建立对应关系，即建立粒子到网格的映射。通过映射得到的相关信息，可以知道一个粒子在哪个网格中，也可以知道在一个网格中有哪些粒子。

上述操作都在本发明所涉及的并行计算系统中的CPU内实现，在得到上述操作的结果后，需要由CPU将相应的信息分配到计算节点内的GPU上。由于GPU间具有并行性，而GPU内的线程间也具有并行性，因此，为了提高计算效率，需要将各个子区域内的信息赋予对应的GPU，将子区域中各个网格内的粒子信息分别赋予对应的线程块(block)。所述的线程块是一个GPU内网格计算的对应单位，即如图3所示，一个线程块对应一个网格，线程块中的一个线程(thread)对应网格中的一个粒子。每个线程块中所开辟的线程数与网格内的粒子数相关，在原则上需要保证一个粒子有一个专门的处理线程。由于粒子在运动过程中会发生位置变化，因此，在一个网格中的粒子数目是动态变化的，所以线程块中所开辟的线程数应具有冗余量，超过网格内粒子数的最大可能值。在前面的说明中提到，根据网格大小和粒子密度，一个网格内的粒子数目在100左右，因此，在GPU线程块内可以开辟160个左右的线程。CPU将信息如前所述分配到GPU内后，还需要为每一个GPU分别分配一个进程做相应的管理。

在为GPU的线程块分配对应网格的粒子后，还需要将粒子的作用势参数、模拟设置信息等信息保存到GPU的显存中。

在背景技术中已经提到，在GPU中采用选择分支的方式对不同作用势进行计算会降低多组分系统离散模拟的计算效率与可靠性。考虑到在一种模拟中，各种粒子间的不同作用势的基本形式相同，有差异的仅是作用势参数。因此，本发明可以将各个作用势参数存储在GPU的一个参数数组中，然后在后续的计算粒子对相互作用时，根据相关粒子的种类属性，选择相应的作用势参数，从而实现粒子间的作用计算。以本实施例所要模拟的三组分系统为例，该系统具有气、液、固三相，因此在气(粒子种类号fos＝1)-液(粒子种类号fos＝2)-固(粒子种类号fos＝3)三相体系中，分子间两两作用共有6种Lennard-Jones势参数。所述势参数具体包括：作用势1，气体-气体间，r_cgg＝2，σ_gg＝0.5，ε_gg＝1，a_gg＝1，b_gg＝0；作用势2，液体-液体间，r_cll＝5，σ_ll＝1，ε_ll＝1，a_ll＝1，b_ll＝1；作用势3，固体-固体间，r_css＝5，σ_ss＝1，ε_ss＝1，a_ss＝1，b_ss＝1；作用势4，气体-液体间，r_cgl＝3，σ_gl＝1，ε_gl＝1，a_gl＝1，b_gl＝0；作用势5，气体-固体间，r_cgs＝3，σ_gs＝1，ε_gs＝1，a_gs＝1，b_gs＝0；作用势6，液体-固体间，r_cls＝5，σ_ls＝1，ε_ls＝1，a_ls＝1，b_ls＝0。其中的r代表作用截断距离。

要保存到GPU显存中的模拟设置信息则包括：质量力g＝2×10^-3，温控采用壁面分子速度标定方法，Y、Z方向为周期性边界，固体粒子由倔强系数C＝75的弹簧“锚”在平衡位置附近运动。

GPU显存除了要存储上述信息外，还要存储前述的粒子到网格映射信息。上述多种信息在GPU中所存储的位置有所不同，其中的粒子到网格映射信息传输到全局内存中，作用势参数以及模拟设置信息存储在常量内存中，利用常量内存的缓存机制来提高作用势计算中对相关势参数的调用速度，提高计算性能。

通过以上操作，为各个GPU分配任务并赋予相应的参数后，就可以在各个GPU的线程块中为对应网格内的粒子执行相应的计算。在计算过程中，由于所要模拟的三组分系统中的各个粒子是随着时间而不断运动着的，因此需要对粒子的位置和速度信息做更新，了解粒子当前所在的位置和速度的大小。在更新粒子位置和速度时，可采用前述参考文献1中所提到的蛙跳方式，即：

$v(t+\mathrm{\Δt}/2)=v\left(t\right)+\frac{f\left(t\right)}{2m}\mathrm{\Δt}$

r(t+Δt)＝r(t)+Δtv(t+Δt/2)

$v(t+\mathrm{\Δt})=v(t+\mathrm{\Δt}/2)+\frac{f(t+\mathrm{\Δt})}{2m}\mathrm{\Δt}$

具体而言，在更新粒子位置和速度时，首先将分子在t时的速度v(t)根据当时的受力f(t)更新到v(t+Δt/2)，从而将其位置r(t)更新到r(t+Δt)；然后按t+Δt时的各分子位置计算受力f(t+Δt)，把速度由v(t+Δt/2)更新到v(t+Δt)，从而完成分子位置、速度的同步。

在更新粒子的位置后，需要重新实现粒子到网格的映射。在前文中已经提到，粒子的位置在运动过程中会发生变化，可能会从一个网格到了另一个网格，因此，需要对粒子到网格的映射关系加以更新。在更新过程中，需要对网格发生变化的粒子重新建立粒子到网格的映射，并删除原网格内的相应映射。在并行计算系统内，各个网格内的粒子映射更新可以在各个线程块内并行进行，即每个网格可同时完成删除已不属于本网格的粒子，并从周边邻居网格接收属于自己的粒子。以一个网格为例，对GPU内的对应线程块如何完成粒子到网格映射的更新进行说明：

步骤a、根据粒子的新位置计算移动后粒子所在网格的坐标；

步骤b、每个网格剔除不属于自己的粒子。本步骤中，首先要将原有的粒子到网格的映射信息进行备份，然后为了防止内存冲突，在线程块内用唯一的活动线程遍历整个网格内的各个粒子的归属情况，若某个粒子所归属的网格的坐标与所在网格不符，则将该网格粒子映射中最后一个实粒子转移到待删除粒子所在位置，而将原最后实粒子位置替代为虚粒子，完成粒子删除操作。

步骤c、网格从周边邻居网格中找出属于自己的粒子，并将这些粒子放入本网格的粒子映射信息中。本步骤中，为了防止内存冲突，在每个线程块中采用唯一的活动线程从步骤b所备份的原有的粒子到网格的映射信息中遍历该网格周边所有邻居网格的所有实粒子新归属，若粒子的新归属坐标与该网格一致，则将粒子插入到该网格的粒子映射信息中。

通过上述的操作，完成了粒子到网格映射的更新，粒子到网格映射的更新主要解决了在同一个子区域内的粒子在不同网格间迁移的问题，但在实际应用中，粒子也有可能从一个子区域移动到另一个子区域，即由一个GPU转到另一个GPU做相关的粒子计算，因此还需要实现在不同GPU间粒子的迁移操作。由于GPU与GPU之间不能直接做数据交换，需要通过CPU的进程完成GPU间数据交换的操作，因此，GPU间粒子迁移的操作也被称为进程间的粒子迁移。为了方便理解，以图4为例，对进程间粒子迁移的操作进行说明。在图4中，根据前述任务划分的需求，整个模拟区域[GlobalLeftBottom，GlobalRightTop]被分割成Nx×Ny个子区域，一个进程Pi，j(进程编号为i+j×Nx)负责一个子区域内([LocalLeftBottom，LocalRightTop])的计算。在前面已经提到，子区域可以进一步划分为网格，在图4中，用小方块表示网格。子区域按照由外到内的顺序，还可以分为外边界区、内边界区以及子区域内部，其中，一个子区域的外边界区是相邻子区域的内边界区，类似的，一个子区域的内边界区是相邻子区域的外边界区。例如，图4中进程P_i，j所对应子区域左侧外边界区是其左侧相邻的进程P_i-1，j所对应的子区域的右侧内边界区。如果某个粒子从进程P_i，j所对应子区域左侧内边界区运动到了左侧外边界区，实际上也就运动到了进程P_i-1，j所对应的子区域的右侧内边界区，因此对该粒子的处理由进程P_i，j所管理的GPU转移到了由进程P_i-1，j所管理的GPU内。

此外，在计算粒子间相互作用力时，由于子区域的边界区域内的粒子需要其他子区域内的粒子信息才能完成完整的受力计算，因此一个GPU还需要得到邻居GPU的边界粒子信息，这种现象也被称为边界粒子传递。边界粒子传递同样需要借助CPU内的进程实现。

上述的进程间粒子迁移与边界粒子传递都可以采用相同的方法实现，因此下面做统一的说明：

步骤A、GPU线程块中的各个线程各自将内边界区域和外边界区域中粒子到网格的映射信息保存到GPU显存的缓冲区内。

步骤B、将显存缓冲区内的信息复制到主存的发送缓冲区内，供进程间通讯传输使用。

步骤C、由CPU实现进程与邻居进程之间数据的发送接收传输。将发送缓冲区内的信息传输到目标进程，并将其他进程传入的信息放置于接收缓冲区内，与所述发送和接收相关的进程均为本进程的邻居进程。

步骤D、将接收缓冲区内信息复制到显存缓冲区内，供GPU下一步操作调用。

步骤E、GPU解析显存缓冲区内的粒子信息，将跨界粒子添加到本进程内边界网格粒子映射中，或者替换外边界网格的粒子映射信息。在跨界粒子添加到对应网格的粒子映射时，只有一个线程在进行写入操作避免内存冲突和数据不完整。

在上述的步骤C中，为了提高通信效率，可以采用现有技术中的shift模式实现相邻进程间的通信。以图4中的空间二维分割为例，shift模式将通信过程分为两个阶段，分别是x方向的传递和y方向的传递，如图5所示。对于空间二维分割而言，shift模式中一个进程只需与相邻四个进程进行通信，与点对点传输模式相比，shift模式通信次数大大减少。

借助CPU完成关于进程间粒子迁移和边界粒子传递的相关操作后，就可以对粒子间作用力进行计算。在计算某一粒子与周边邻居粒子间的作用力时，该粒子的信息需要被频繁访问，因此将粒子信息(前述的粒子到网格映射信息中已经包含了粒子信息)放置在GPU的共享内存中以提高访问速度进而提高计算性能。对粒子间作用力的计算通常对一个线程块内的所有粒子同时进行，其实现步骤如下：

步骤1、GPU在执行的内核(kernel)中为每个线程块开辟共享内存，以存放该线程块对应的网格以及它的邻居网格，所存放的网格信息中包括粒子信息；

步骤2、GPU从全局内存中读入网格信息后复制到对应的共享内存；

步骤3、在共享内存内，计算线程块所对应网格内的粒子与一个邻居网格内粒子的所有作用力，在计算过程中邻居网格内的某一粒子能被该线程块内的所有线程同时访问；

步骤4、在计算完一个线程块所对应网格内的粒子与一个邻居网格内的粒子间的相互作用力后轮换邻居网格，重复上述的步骤3，直到执行完与周边所有邻居粒子(几百～几千个)间的作用力计算后才将线程块对应网格内的粒子的信息写回全局内存。

以上就是粒子间作用力计算的基本实现步骤。上述步骤在实现粒子与其他粒子间的作用力时，参照背景技术中所提到的相关内容可以知道，在计算粒子间作用力时需要用到作用势，而计算作用势时又需要根据粒子的类型选择适当的作用势参数。在前面的说明中，已经对本实施例所模拟的三组分系统中所可能包含的势参数做了列举，在此则决定具体选用哪个势参数。在如下的表1中，对本实施例中三组分系统的势函数参数选择方式做了说明：

表1

根据上表中的内容，结合下列公式，可以知道某线程计算两个粒子间的作用力时所选用的势参数：

PotentialType＝(part1.fos＝＝part2.fos)*part1.fos+(part1.fos！＝part2.fos)*(part1.fos+part2.fos+1)

在得到势参数后，就可以对粒子间的作用力进行计算。具体的计算过程在背景技术中已经做了较为详细的说明，此处不再重复。

在得到粒子所受到的作用力后，可以根据作用力再次更新粒子的速度，对粒子速度更新的实现可参照前面的说明。然后对三组分系统的模拟是否结束进行判断，若没有结束，则重新执行上述从更新粒子位移起的各个操作步骤，直到模拟结束后将相关信息从并行计算系统输出。

以上是对在带有GPU的并行计算系统上实现三组分系统模拟的基本流程。采用上述方法可以取得良好的模拟效果。对上述实施例中所模拟的三组分系统的模拟结果在图6中做了展示。其中，图中的体积较大的絮状物(在图6(a)中为体积较大的圆球)表示气泡，体积较小的圆球状物表示固定颗粒，填充在空间内的其余部分为液体。在图6(b)中t＝300时，气泡第一次接触固体颗粒，在图6(c)中t＝500时，固体颗粒穿过气泡；而在图6(d)中t＝900时，由于固体颗粒周围流场非对称，气泡受力偏倚，之后沿一侧随液体流动。在图6(e)中t＝1250时，气泡再次与固体颗粒接触，由于气液界面张力的限制而未被固体颗粒穿透，而后气泡贴在壁面运动。此现象与人们设想的矿物浮选等过程中气固接触的一些典型过程相符，因此采用上述方法所实现的模拟结果较为真实可信。

经过测算，本实施例中的GPU并行计算系统发挥的计算性能约为450Gflops，是单个CPU可发挥性能的200倍左右，可见本发明的方法在计算效率上也有很大的提高。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (9)

1.一种基于图形处理器的多组分系统离散模拟计算方法，其特征在于，该方法在由多个计算节点组成的并行计算系统上实现，所述计算节点上安装有图形处理器GPU；该方法包括：

步骤1)、对所要模拟的多组分系统做空间分解，在所述空间分解中，将所述多组分系统的空间分为多个子区域，并将所述子区域分为网格；

步骤2)、建立所述多组分系统中的粒子与所述步骤1)所得到的网格之间的映射关系；

步骤3)、将所述多组分系统中的粒子信息按照所述空间分解的结果分别分配给所述多个计算节点上的GPU内的线程块，并将粒子的作用势参数、模拟设置信息保存到所述GPU的内存中；

步骤4)、在所述GPU上计算粒子所受到的作用力，并根据所述作用力对粒子的速度与位置进行更新；

步骤5)、重复对粒子的上述计算过程，直到对所述多组分系统模拟结束后将模拟结果从所述并行计算系统中输出。

2.根据权利要求1所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，所述的步骤4)包括：

步骤4-1)、在所述GPU的各个线程块上并行计算各个粒子的当前速度与当前位置，根据计算结果实现对粒子位移以及粒子速度的更新；

步骤4-2)、根据步骤4-1)的更新结果在各个GPU内实现所述粒子到所述网格映射的更新；

步骤4-3)、利用所述计算节点上的CPU进程实现GPU间的信息传递，实现粒子在不同GPU间的迁移以及边界粒子传递；

步骤4-4)、根据所述粒子的类型，选择相应的粒子作用势参数，实现粒子间作用力的计算；

步骤4-5)、根据所述粒子间作用力的计算结果，再次对所述粒子的速度做更新操作。

3.根据权利要求1或2所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，在所述的步骤1)中，按照所述并行计算系统中的GPU个数实现对所述子区域的划分，使得所述子区域的个数与所述GPU的个数相等；在所述子区域内按照作用截断距离信息划分网格。

4.根据权利要求1或2所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，在所述的步骤3)中，将所述多组分系统中的粒子信息按照所述空间分解的结果分别分配给所述多个计算节点上的GPU内的线程块时，将所述多组分系统的一个子区域分配给一个所述的GPU，将所述子区域内的一个网格分配给所述GPU内的一个线程块，所述网格内的一个粒子由所述线程块内的一个线程做专门的处理。

5.根据权利要求1或2所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，在所述的步骤3)中，所述的作用势参数包括作用截断距离r、作用特征长度σ、作用势能ε、用于表示斥力项相对大小的系数a、用于表示引力项相对大小的系数b在内的多个参数。

6.根据权利要求2所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，所述的步骤4-2)包括：

步骤4-2-1)、根据粒子的新位置计算粒子在运动后所在网格的坐标；

步骤4-2-2)、各个网格备份原有的粒子到网格的映射信息，并从所述映射信息中找出网格坐标与当前网格坐标不相符合的粒子，删除该粒子；

步骤4-2-3)、网格根据周边邻居网格所备份的粒子到网格的映射信息，从周边邻居中找出属于自己的粒子，并将这些粒子放入本网格的粒子映射信息中。

7.根据权利要求2所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，所述的步骤4-3)包括：

步骤4-3-1)、GPU线程块中的各个线程分别将位于内边界区域和外边界区域中的各个粒子的到网格的映射信息保存到所述GPU内存的缓冲区内；

步骤4-3-2)、将GPU内存缓冲区内的信息复制到对所述GPU进行管理的CPU进程所对应的发送缓冲区内；

步骤4-3-3)、由所述CPU进程将所述发送缓冲区内的信息传输到目标邻居进程的接收缓冲区，将其他邻居进程传入的信息放置于本进程的接收缓冲区内；

步骤4-3-4)、将所述接收缓冲区内的信息复制到所述GPU的内存缓冲区；

步骤4-3-5)、GPU解析内存缓冲区内的粒子信息，将跨界粒子添加到内边界网格粒子的映射信息中，或者对位于外边界网格的粒子的映射信息做替换。

8.根据权利要求7所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，在所述的步骤4-3-3)中，CPU进程与邻居进程间的通信采用shift模式。

9.根据权利要求2所述的基于图形处理器的多组分系统离散模拟计算方法，其特征在于，所述的步骤4-4)包括：

步骤4-4-1)、GPU在内核中为一个线程块开辟共享内存；

步骤4-4-2)、在所述的共享内存内存放该线程块所对应的网格的信息以及该网格的邻居网格的信息；

步骤4-4-3)、在共享内存内，根据粒子类型，选择相应的粒子作用势参数，计算线程块所对应网格内的粒子与一个邻居网格内粒子的所有作用力；

步骤4-4-4)、在计算完一个线程块所对应网格内的粒子与一个邻居网格内的粒子间的相互作用力后轮换邻居网格，重复上述的步骤4-4-3)，直到执行完与周边所有邻居粒子间的作用力计算后才将线程块对应网格内的包括作用力大小的粒子信息写回。

Images