CN102760116A - 一种基于硬球模型的并行计算方法 - Google Patents

Abstract

本发明提供一种基于硬球模型大规模并行计算的方法，该方法将所要模拟的系统划分为多个子区域，同时将所述子区域划分为多个内部区元胞网格和边界区元胞网格，并且在所述子区域外层建立子区域的虚拟区元胞网格；所述内部区元胞网格对应事件驱动区域，所述边界区元胞网格和虚拟区元胞网格对应时间驱动区域；在给定的时间驱动时间步中，更新所述事件驱动区域和时间驱动区域内粒子的状态。本发明克服了硬球模型事件驱动算法难以并行且现有并行算法由于反复回调而导致的低效困难，同时也解决了时驱算法由于计算量大而造成的计算效率低下问题，在计算效率和并行扩展性上有了很大的提高。

Description

一种基于硬球模型的并行计算方法

技术领域

本发明涉及离散粒子模拟并行计算，特别涉及到硬球模型的大规模可扩展并行计算。

背景技术

在颗粒流体系统中，普遍存在着流动、传递、反应等多种过程的耦合，如化工、能源等领域中涉及颗粒和流体的流态化系统，使得系统呈现出显著的非平衡性和多尺度结构等。目前对这些系统的研究虽然已经有一些成果，但总体而言，要严格描述这些系统的多尺度结构还是十分困难的，究其原因还是其机理太过复杂，从而难以进行深层次的研究。

当前，对于颗粒流体系统等领域的研究主要还用采取连续介质模型，通过欧拉(Euler)坐标系进行直接的数值模拟，其尺度可追溯到连续介质模型所能成立的最小尺度。但是此类方法在处理介观结构和稀薄气体等方面存在很大困难，此时连续介质假设已经不再成立，需要采用离散的方法来探索和研究。从本质上讲，自然界中的物质均是离散的，因此采用离散的方法对系统进行研究也就是很自然的事情。

离散模拟(Discrete simulation)的研究方法采用拉格朗日(Lagrange)坐标系，着眼于构成物质体系的质点，将系统离散为经典力学下的离散质点，直接描述系统的底层物理结构和过程。此类方法不依赖于宏观的纳维叶-斯托克斯(Navier-Stokes)方程等唯象描述，深入而可靠地反映了系统的基本机理，被越来越多地应用到两相体系的模拟研究中。

一般来说，处理离散粒子间的相互作用及其变化的方法主要有两大类：软球(Soft Sphere)模型和硬球(Hard Sphere)模型(二维模拟中则称之为软碟Soft Disk和硬碟Hard Disk)。软球模型通过弹性、阻尼以及滑移的力学机制考虑颗粒碰撞时的相互作用，力的计算和系统的推演，具有很好的并行性。而硬球模型则认为粒子间的作用是瞬时的二元碰撞，碰撞后粒子的运动速度按经典力学计算，特别适合气体尤其是稀薄气体流动的模拟计算，从而得到了广泛的应用。

由于研究体系的不断增大，模拟的粒子数目也越来越多。为了解决应对巨大的计算量，需要进行并行计算。不过由于硬球模型采用的是事件驱动算法(Event-driven algorithm)，虽然单机的计算效率相当优越，但是其内在性质使之难以实现真正意义上的大规模并行计算。目前虽然存在硬球模型事驱算法的并行方案(参见参考文献1：Marín M.，Computer PhysicsCommunications，102(1-3)，81-96，1997；参考文献2：Miller S.，Luding S.，Journal of Computational Physics，193(1)，306-316，2004)，并初步解决了此问题。但是当研究的规模变大以后，算法中的因果错误出现的几率会越来越大，程序也将会频繁地进行回调，从而降低了算法的效率，并且使得算法并不能真正意义上实现大规模可扩展并行计算。

不过基于硬球模型的拟颗粒模型(Pseudo-particle Modeling，PPM)(参见参考文献3：Ge W.，Li J.H.Pseudo-Particle Approach to Hydrodynamics ofParticle-Fluid Systems.In：Proceedings of the Fifth International Conferenceon Circulating Fluidized Beds，Beijing，Science Press，1996)，则不存在并行的问题。拟颗粒模型将气体处理为大量的硬球粒子，按时间驱动的方式推演。拟颗粒方法结合了分子动力学中硬球模型模拟和直接蒙特卡洛模拟两种模拟方法的优点，因此具有较好的微观表现力和计算效率，且基于时间驱动的本质，使得拟颗粒模型具有良好的并行性。但是，由于拟颗粒模型在模拟过程中每步都要判断粒子之间碰撞关系，对于气体尤其是稀薄气体这种碰撞概率比较低的体系则浪费了大量的计算机时，从而造成模拟的效率过于低下。

发明内容

本发明的目的是克服现有事件驱动算法难以并行以及现有硬球时间驱动算法计算效率低下的不足，针对采用硬球模型的离散粒子系统，提供一种基于时间驱动算法和事件驱动算法耦合模式的大规模可扩展并行算法。

为了实现上述目的，本发明提供了一种基于硬球模型的并行计算方法，该方法包括以下步骤：

步骤1)将所要模拟的粒子系统信息读入并行计算系统中；

步骤2)将所要模拟的系统划分为多个子区域，同时将所述子区域划分为多个元胞网格并且在所述子区域外层建立子区域的虚拟区元胞网格，所述子区域元胞网格分为边界区元胞网格和内部区元胞网格；其中，所述内部区元胞网格对应事件驱动区域，所述边界区元胞网格和虚拟区元胞网格对应时间驱动区域；

步骤3)建立所要模拟的系统中粒子与元胞网格的映射关系；

步骤4)在给定的时间驱动时间步中，更新所述事件驱动区域和时间驱动区域内粒子的状态。

其中，可以重复执行步骤4)直至系统达到稳定状态或者满足用户要求，然后将模拟结果从并行计算系统中输出。

根据本发明实施例的并行计算方法，在所述步骤2)中每个子区域分别对应所述并行计算系统中的一个进程；所述虚拟区元胞网格和边界区元胞网格的层数均不少于2层，并且虚拟区元胞网格的层数不超过边界区元胞网格的层数。

根据本发明实施例的并行计算方法，在所述步骤3)中建立所要模拟系统中粒子与元胞网格的映射关系是指根据粒子的位置将粒子以链表方式分布至步骤1)中所述子区域元胞网格中。

根据本发明实施例的并行计算方法，在所述步骤4)中所述给定的时间驱动时间步Δt的大小应满足：

|v₁-v₂|Δt≤min(d/2，l₁₂-r₁-r₂)

其中，v₁、v₂是两个粒子的速度，l₁₂是粒子之间的距离，d是粒子的平均直径，r₁和r₂是粒子的半径。

根据本发明实施例的并行计算方法，在所述步骤4)中更新所述事件驱动区域内粒子的状态包括以下步骤：

步骤4-1)建立所述子区域内部区元胞网格的事件网格；

步骤4-2)将所述子区域中每个事件网格中的粒子的事件按发生的先后顺序存储在事件网格的一个链表中；

步骤4-3)找到所述子区域内部区元胞网格所对应的所有事件网格中最优先发生事件；

步骤4-4)如果所述最优先发生事件为穿越元胞网格事件则更新其时间和位置，并更新粒子所处元胞网格所在的事件网格的事件链表；同时分析粒子新位置所处区域，如果粒子处于时驱区域则更新粒子时间至时驱时间步并继续更新粒子的位置，如粒子的最新位置在事驱区域，则更新该最新位置所在元胞网格的事件网格的事件链表；

步骤4-5)如果所述最优先发生事件为碰撞事件则更新碰撞粒子对的时间、位置和速度；分析碰撞粒子对新位置所处元胞网格，并更新粒子链表发生更改的元胞网格所处事件网格中的事件链表；

步骤4-6)重复执行步骤4-3至4-5，直到处理完给定时间步内事驱区域所有事件。

根据本发明实施例的并行计算方法，在所述步骤4)中更新所述时间驱动区域内粒子的状态包括以下步骤：

步骤4-7)更新所述边界区元胞网格中的粒子的位置信息，并完成粒子的进出进程操作；

步骤4-8)、处理所述边界区元胞网格中粒子之间的碰撞，以及边界区元胞网格中粒子和内部区元胞网格中粒子的碰撞，并更新该粒子和被碰撞粒子的速度；

步骤4-9)、处理本子区域边界区元胞网格中粒子与相邻子区域边界区元胞网格中粒子之间的碰撞，并更新碰撞粒子对的速度。

根据本发明实施例的并行计算方法，在所述步骤4-7)中更新粒子的位置信息包括以下步骤：

步骤4-7-1)、如果粒子的位置仍然处于原先网格，则粒子所属网格信息不变；

步骤4-7-2)、如果粒子的位置不在原先网格，但粒子的新位置仍处于所述边界区和虚拟区元胞网格中，则将粒子移至新位置所处网格中；如果粒子的位置不在原先网格，且粒子的新位置处于内部区元胞网格中，则将粒子移至新位置所处网格中，并更新新位置所处网格所在的事件网格的事件列表；

步骤4-7-3)、将虚拟区元胞网格中的粒子通信至相邻子区域，并从虚拟区元胞网格中删除粒子的信息。

根据本发明实施例的并行计算方法，所述步骤4-9)是依次按以下顺序处理处理本子区域边界区元胞网格中粒子与位于其它子区域边界区元胞网格中粒子之间的碰撞的：首先处理具有最少相邻子区域的边界区元胞网格，接着处理具有次少相邻子区域的边界区元胞网格，直到处理完具有最多相邻子区域的边界区元胞网格为止。

根据本发明实施例的并行计算方法，所述步骤4-9)处理本子区域边界区元胞网格中粒子与相邻子区域边界区元胞网格中粒子之间的碰撞包括以下步骤：

a.将本子区域与高位子区域相邻的边界区元胞网格粒子的信息传递至所述高位子区域的虚拟区元胞网格中；

b.本子区域对应的进程计算与低位子区域相邻的边界区元胞网格中粒子和虚拟区元胞网格中粒子之间的碰撞，并更新碰撞粒子的状态；

c.将已碰撞且更新状态后的虚拟区的粒子信息反馈至低位子区域边界区对应的元胞网格之中，并从虚拟区中删除粒子的信息；

其中，所述高位子区域与低位子区域是按照步骤2)中所划分的子区域的位置而确定的。

与现有技术相比，本发明的优点在于：耦合了时间驱动算法和事件驱动算法各自的优点，利用时间驱动算法易于并行的特点，克服了事件驱动算法固有的难以并行的缺陷，解决了时驱算法计算效率低下的问题。并利用并行机群系统对大量数据进行并行处理，加快了整体的运算效率，扩大了计算规模。另外，还提出用于硬球模型高维系统高维分割情形下的通信方式，解决了一般通信方式中由于硬球碰撞时序不一致、而导致系统能量在碰撞处理过程中不守恒的情形。

附图说明

图1为根据本发明实施例的区域分解示意图；

图2为根据本发明实施例的子区域内部区和边界区元胞网格划分以及扩展的虚拟区网格划分示意图；

图3为根据本发明实施例的粒子在元胞网格中的存储以及粒子穿越网格示意图；

图4为根据本发明实施例的对于搜索某一进程的事件驱动区域中最优先发生事件及其更新的二叉树法示意图；

图5为在邻居进程间所采用的Shift通信模式的一个实施例的示意图；

图6为根据本发明实施例的二维模拟系统二维分割边界区域第一次通信示意图；

图7为根据本发明实施例的二维模拟系统二维分割边界区域第二次通信示意图；

图8为根据本发明实施例的三维模拟系统三维分割边界区域的通信示意图；

图9是描述基于硬球模型的大规模并行计算方法的一个实施例的流程图；

图10为根据本发明实施例的基于硬球模型的大规模并行计算方法的计算效率比较图；

图11为根据本发明实施例的基于硬球模型的大规模并行计算方法的计算过程比较图。

具体实施方式

为了使本发明的目的，技术方案及优点更加清楚明白，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在本发明的实施例中，在模拟稀薄气体粒子系统时，以CPU为计算部件的计算机系统为基础实现模拟过程。模拟中，粒子的模型为硬球模型。虽然硬球模型事驱算法单机的效率已经很高，但是为了更好地对大量粒子的模拟，在本发明的实施例中，采用并行计算系统实现对大规模硬球粒子系统的模拟。所述的并行计算系统由多个计算节点组成，在每个计算节点上都采用多个多核处理器，不同计算节点间通过高速网络连接。在本发明的一个实施例中，系统的计算节点具有如下配置：每个节点有两个型号为Xeon 5430的CPU，每个CPU有四个计算核心，每个节点有容量为16GB大小的DDR2-667内存，不同节点间通过InfiniBand网络连接起来。在本实例中，采用了2个上述节点，每个节点使用了8个进程，共计16个进程。本领域的普通技术人员应该理解此处并行系统的配置仅起示意说明的作用，在其他实施例中可以采用现有的各种并行计算系统的硬件配置。

在本发明的实施例中所涉及的硬球作用模型和事件驱动算法是现有技术，为了方便理解以及下文中说明的需要，首先对硬球作用模型和事件驱动算法加以说明。

硬球模型是计算两个粒子之间的作用势的模型，其表达式如下：

$U\left(l_{\mathrm{ij}}\right)=\left\{\begin{array}{c}\∞,l_{\mathrm{ij}}\≤r_i+r_j\\ 0,l_{\mathrm{ij}}>r_i+r_j\end{array}\right.$

上述公式中，U表示粒子之间的势能，l是粒子之间的距离，r是粒子的半径，i，j表示粒子的标号。

关于硬球模型的具体内容在参考文献4“Alder B J，Wainwright T E.The Journal of Chemical Physics，27(5)，1208-1209，1957”中有详细说明。

事件驱动算法中，两个粒子碰撞时刻为：

$t_c=t_1+\frac{1}{v_{12}^2}[-P_{12}\·v_{12}-\sqrt{P_{12}\·v_{12}-(P_{12}^2-{(r_1+r_2)}^2)v_{12}^2}]$

上述公式中，t₁是1号粒子所处时刻，v₁₂是两个粒子的速度差，P₁₂是两个粒子的位移差，r是粒子的半径，下标1、2是粒子的标号。

实例中，两个粒子碰撞后的速度为：

$v_{\mathrm{1,2}}^{\′}=v_{\mathrm{1,2}}m\frac{(1+e)m_{\mathrm{2,1}}}{m_1+m_2}\frac{(v_1-v_2)\·(P_1-P_2)}{{|P_1-P_2|}^2}(P_1-P_2)$

上述公式中，e是恢复系数，v’是粒子碰撞后的速度，v是粒子碰撞前的速度，m是粒子的质量，P是粒子的位置，下标1、2是粒子的标号。

关于硬球事件驱动算法，参考文献1中有详细说明；关于硬球时驱算法拟颗粒模型，参考文献3中有详细说明。

在本发明的实施例中，提供了一种基于时间驱动算法和事件驱动算法耦合模式的大规模可扩展并行方法，在上文所述的由多个计算节点组成的并行计算系统上实现，主要包括下面的步骤：

1、模拟区域的空间分解及并行任务的划分

将所要模拟的粒子系统信息读入到前述的并行计算系统中，形成模拟区域。同时，在所述的计算机系统中实现对该模拟区域的空间分解，以完成并行任务的划分。图1为根据本发明实施例的区域分解示意图，其中，对所要模拟的粒子系统作区域空间分解，将所述粒子系统的空间划分为多个子区域，每个子区域应该对应于一个计算进程。在一些实施例中，所划分子区域的个数通常就是并行计算系统的计算核心的个数，每个计算核心对应一个子区域构成一个计算进程。在本实施例中，按照上述并行计算系统中的进程数可将模拟区域划分为16个子区域。在其他实施例中，可以按用户的需求或当前系统实际可利用资源的情况来划分子区域。仍在其他实施例中，通过输入模拟参数来指定所划分的子区域的个数。

接着，采用元胞列表法(参考文献4：Hockney R.W.，Goel S.P.，Journalof Computational Physics，14(2)，148-158，1974)，将每个子区域划分为许多尺寸相等的元胞网格。在一些实施例中，元胞网格的尺寸大小l_m可以是满足事驱算法最高效的网格尺寸，根据文献5：Mauricio Marin.Journal ofComputational Physics，109，306-317，1993的图4.1回归可得，该尺寸近似满足：

$l_m=\left\{\begin{array}{cc}l/\sqrt{2},& l/d\&GreaterEqual;\sqrt{2}\\ d,& l/d<\sqrt{2}\end{array}\right.$

其中，l是所述系统中粒子间平均距离，d是粒子的直径。在本实施例中，为了更好地划分网格和方便说明，粒子的无量纲质量和半径均为1，填充率为0.0873，则l约为6.0。l_m取为和l相等。此时算法的效率会有所降低，但由于气体比较稀薄，因此降低不是很明显。在其他实施例中，元胞网格的尺寸也可以根据用户指定的参数进行设置。

不同于参考文献4的单节点的元胞列表法，在本发明的实施例中由于采用了多节点来进行并行计算，所以将子区域元胞网格区分为内部区元胞网格和边界区元胞网格；并且在每个子区域的边界区元胞网格外层，建立子区域的虚拟区元胞网格。

图2为根据本发明实施例的子区域内部区和边界区元胞网格划分以及扩展的虚拟区网格划分示意图。其中虚拟区元胞网格和边界区元胞网格的层数均不应少于2层。在一些实施例中，虚拟区元胞网格的层数等于边界区元胞网格的层数，在其他实施例中，虚拟区元胞网格的层数小于边界区元胞网格的层数。仍在其他实施例中，边界区元胞网格和虚拟区元胞网格的层数可以由用户指定，但均不应小于2层。在一个实施例中，虚拟区元胞网格尺寸与在所述子区域内划分的元胞网格相同。在另一实施例中，虚拟区元胞网格的尺寸可以由用户指定。从上文描述可以看出，虚拟区元胞网格实际上对应于邻近子区域的边界区元胞网格或相邻子区域的边界区元胞网格的一部分。

进一步参考图2和图5(a)，本实施例中所模拟的粒子系统为二维系统，划分子区域后，每个子区域有一个或多个相邻子区域。在二维分割的模拟系统中，每个子区域可以有两个或四个相邻子区域。在横向(x方向)上，可以将在本子区域的左边的相邻子区域称为是低位子区域，在本子区域右边的相邻的子区域称为是高位子区域。类似地，在竖向(y方向)上可以将位于本子区域的下面的相邻子区域称为是低位子区域，位于本子区域的上面的相邻的子区域称为是高位子区域。位于X方向起点和Y方向起点相交的位置的子区域可以由用户指定或由系统缺省设置。由于每个子区域对应并行计算系统中的一个计算进程。所以根据所划分的子区域的位置对子区域所对应的计算进程指定进程标号P_i，j，P_i-1，j为P_i，j的x低位进程；P_i+1，j为P_i，j的x高位进程；P_i，j-1为P_i，j的y低位进程P_i，j+1为P_i，j的y高位进程。进程之间的通信次序应满足本申请中下面所讨论的进程通信次序。

在对模拟区域作区域空间分解以后，还需要建立粒子信息和元胞网格之间的映射。

图3为根据本发明实施例的粒子在元胞网格中的存储以及粒子穿越网格示意图。其中，并行计算机系统根据模拟区域中的粒子所在的位置，将粒子系统中粒子信息按照所述空间分解的结果分别保存到所述多个计算机进程的系统内存中，同时根据粒子的位置将粒子布置到所在位置的元胞网格之中，建立粒子和元胞网格的映射。粒子在元胞网格中通过链表方式存储，元胞网格中的每个粒子被保存为链表中的一项，当粒子离开该元胞网格时，从链表中删除相应的项；当粒子进入一个元胞网格时，在该元胞网格的粒子的链表中加入相应的项。这样通过映射关系可以知道一个粒子在哪个元胞网格之中，也可以知道一个元胞网格中有哪些粒子，从而将粒子分割在不同的子区域的元胞网格之中。

2、耦合时驱和事驱算法的并行计算方法

在本发明的实施例中采用了耦合时间驱动和事件驱动的并行计算方法，对不同位置的粒子采用不同的算法：子区域的内部区元胞网格及其内部粒子对应为事件驱动；子区域的边界区元胞网格和虚拟区元胞网格以及这些元胞网格内部粒子对应为时间驱动。每个子区域与一个计算进程相对应，在多个计算进程上并行地计算粒子之间的碰撞，更新粒子的时间、位置等信息。为了方便，在下面的叙述中有时会将计算进程所处理的对应子区域的元胞网格的简称为本进程的元胞网格，将子区域所对应的计算进程简称为本计算进程等，但本领域的普通技术人员通过上文的描述应该理解一个子区域对应一个计算进程，计算进程所处理的是对应子区域中的粒子。具体计算步骤如下：

步骤1，在给定的时间驱动时间步中，更新位于事件驱动区域中的粒子的状态。

所述给定的时间驱动时间步Δt的大小应满足：

|v₁-v₂|Δt≤min(d/2，l₁₂-r₁-r₂)

其中，v₁、v₂是两个粒子的速度，l₁₂是粒子之间的距离，d是粒子的平均直径，r₁和r₂是粒子的半径。时间步长取得过小，则系统的演化也会变慢，也间接降低了算法的效率；取得太大虽然会增快系统的演化，间接提升算法的效率，但是会导致模拟的结果不是很准确。对于本实施例中的系统，时间步长取为1.0，该值兼顾了系统演化的速度和模拟的准确性。时间步长具体取值方法可参考文献3。

对某一进程对应的子区域的内部区而言，该部分采用的事件驱动算法。参考文献1中对事件的处理以粒子为单位处理，对每个粒子将要发生的事件(与邻近粒子的碰撞或穿越元胞网格)按时间顺序排列(即LMA算法，Local Minima Algorithm，见参考文献1)。与之不同的是，在本发明的实施例中对事件的处理以区域为单位进行，因此可以称为是对区域的LMA算法，即，对子区域内部区中的事件也采用划分网格的方式，将每个进程对应的子区域划分为很多事件网格。这样对事件的处理实际上是以事件网格为单位进行，这样的处理避免了对粒子LMA算法(参考文献1)由于粒子进出进程而导致之后提及的二叉树中算法中节点频繁的删除与添加操作，另一方面也方便了区域分解的并行算法对事件的处理以及时驱事驱耦合算法的实现。为了计算方便，在本实施例中，对于每个进程，所划分的事件网格和上文所述的本地元胞网格重合，即事件网格的大小取为与元胞网格尺寸相等。在其他实施例中，可以根据实际需求划分不同尺寸的事件网格。将进程中每个事件网格中的粒子的事件按发生的先后顺序存储在事件网格的一个链表中；本进程中所有事件网格中最优先发生事件构成该进程的事件列表。

图4为根据本发明实施例的对于搜索某一进程的事件驱动区域中最优先发生事件及其更新的二叉树法示意图；其中，如图4(a)所示，在该进程事件驱动区域内分布的8个事件网格中最优先发生事件的时间依次分别为0.17、0.11、0.23、0.13、0.15、0.22、0.34、0.23，按照二叉树法搜索出事件网格2中最优先发生的时间为0.11的事件作为该进程的最优先发生事件。处理完事件网格2中发生时间为0.11的事件后，更新事件网格2中的事件链表，将事件网格2中的最优先发生时间设置为更新事件链表后的最优先发生事件的时间0.42，再次利用二叉树法更新该进程最优先发生事件，选出具有最小发生时间0.13的事件网格4中的最先发生事件作为当前进程最优先发生事件(如图4(b)所示)。

根据图4所示的方法从本进程的事件列表中选出的本进程的最优先发生事件，然后根据最优先发生事件的类型进行相对应的处理，如果为穿越元胞网格事件则更新粒子的时间和位置；如果碰撞事件则更新碰撞粒子对的时间、位置和速度。

当搜索到的最优先发生事件是碰撞事件时，则将两个粒子的时间提升至即将碰撞的时刻，并且将它们的位置也更新到碰撞时刻所处位置，此后两个粒子发生碰撞，更新两个粒子的速度。分析碰撞粒子对新的位置所处的元胞网格，如果两个粒子在同一元胞网格，则只需更新该元胞网格所在事件网格中的事件链表；如果两个粒子不在同一元胞网格，则必须更新两个粒子各自所处元胞网格所在事件网格中的事件链表。

值得指出的是，事件驱动算法在处理位于事件驱动区域中的某粒子的碰撞中，即使该粒子所处元胞网格与边界区相邻，搜索该粒子和周边粒子的碰撞关系时，也不搜索位于边界区元胞网格中的粒子。事件驱动区域中粒子和时间驱动区域中粒子之间的碰撞在之后的步骤中处理。

当搜索到的最优先发生事件是粒子穿越网格事件时，则首先将粒子的时间和位置更新到事件发生时刻的状态，并更新原先元胞网格所在事件网格中事件的链表，然后分析粒子的新位置所在的元胞网格。如果粒子的新位置所在元胞网格在事驱区域，则需更新粒子新位置所在的元胞网格所在事件网格的事件链表。如果粒子新位置所在元胞网格在时驱区域，则应更新粒子的时间至给定的时间驱动的时间步，并更新此时粒子的位置。如果粒子的最新位置所在元胞网格和之前的新位置所在元胞网格不是同一个网格，则应将粒子移至最新元胞网格之中，并分析该最新元胞网格所在区域。如其所在区域为事件驱动，还需更新该元胞网格所在事件网格的事件列表。

步骤2，在给定的时间驱动时间步中，更新位于时间驱动区域中的粒子的状态。

首先，更新粒子的位置信息，并完成进程间的粒子迁移。如果粒子的位置仍然处于原先网格，则粒子所属网格信息不变；如果粒子的位置不在原先网格，但粒子的新位置仍处于时驱区域，则将粒子移至新位置所处网格中；如果粒子的位置不在原先网格，且粒子的新位置处于事驱区域，则将粒子移至新位置所处网格中，并更新该粒子的新位置所处网格所在的事件网格的事件链表；如果该粒子脱离本进程计算区域，则将其移至虚拟区所对应元胞网格位置。

在进行碰撞处理之前，还需要将上述脱离该进程的粒子通过Shift通信方式(参见参考文献6：Plimpton，Journal of Computational Physics，117(1-19)，1995；参考文献7：Clark T W，Hanxleden R V，McCammon J A andScott L R.Paralleling molecular dynamics using spatial decomposition.Proceedings of the scalable high performance computing conference，Knoxville，TN，1994，95-102)通信至相应的邻近进程，并从该进程中删除相应的粒子的信息。也就是将虚拟区中的这些粒子，通信至相邻进程的对应的本地元胞网格之中，并在本进程虚拟区中删除这些粒子的信息。在又一个实施例中，也可以采用现有的为本领域普通技术人员所熟知的其他并行通信技术来在进程间传递粒子信息。在其他实施例中，也可以使用点对点通信方式，但该通信方式会增加通信的次数，从而造成并行算法效率的下降。但是应指出在上述的实施例中进程间通信的次序应满足本申请下面所讨论的通信次序。

其次，处理时间驱动区域中粒子的碰撞。时间驱动区域包括边界区和虚拟区。此时，对于本计算进程而言，其虚拟区中已不存在粒子，要处理的碰撞主要为本进程内边界区内粒子之间的碰撞、本进程内的边界区中粒子和内部区中粒子之间的碰撞，以及在处理完本进程内的时间驱动区域中粒子的碰撞之后，处理相邻进程之间边界区中粒子的碰撞。而实际上，本进程的虚拟区所对应的元胞网格中也可能会有粒子，因为从上文描述可以看出，本进程的虚拟区元胞网格实际上还对应于邻近进程中的边界区元胞网格，可见对于虚拟区内部粒子之间的碰撞，可通过对各个进程的边界区内部粒子之间的碰撞来处理。不同位置的粒子间碰撞有不同的处理方式。

如果是本计算进程的边界区网格中粒子之间的碰撞，则分析两个粒子之间的距离，如果距离小于等于两个粒子的半径之和且两粒子间有相互靠近的趋势，则两个粒子之间发生碰撞，粒子碰撞后的速度和上文介绍的事件驱动算法中的计算公式相同。

如果是本计算进程的边界区网格中粒子和内部区中粒子的碰撞，则首先将内部区中的粒子的时间提升至时间驱动的时刻，并分析两者之间的距离关系。如果两个粒子之间的距离小于两者半径之和且两粒子有相互靠近的趋势，则两个粒子发生碰撞，碰撞后的速度和上文介绍的事件驱动算法中的计算公式相同。如果两粒子发生碰撞，还必须分析刚提升时间的粒子提升后的位置。如果该粒子移入新的内部区元胞网格，则必须更新该元胞网格所处事件网格中的事件链表，否则只需更新其原先元胞网格所处时间网格中的事件链表。

处理完上述本进程内的边界区的粒子碰撞之后，开始处理相邻进程之间边界区粒子之间的碰撞。

本计算进程中边界区网格中粒子与位于其它计算进程边界区网格中粒子之间的碰撞，由于子区域的边界区域内的粒子需要其它子区域内的粒子信息才能完成的碰撞处理，因此每个进程还需要得到邻居的进程的边界粒子信息才可以完成，这种现象称为边界粒子碰撞处理。但是由于硬球碰撞的瞬时性，硬球之间的碰撞处理顺序应该对于每个进程是一致的，否则会造成系统能量的不守恒，这就是硬球碰撞处理的一致性操作。为了解决这个问题，在本发明的实施例中对于边界区域的粒子的碰撞处理实行分类处理，并采用单向传递原则。下面结合附图具体介绍边界区粒子的通信方式。

图5为在邻居进程间所采用的Shift通信模式的一个实施例的示意图；其中，图5(a)描述了Shift模式的通信的次序。首先接收来自x方向邻近低位进程的粒子信息，处理本进程边界区粒子碰撞，并将更新后的粒子信息反馈给x方向邻近低位进程，类似地处理y方向的进程间的粒子信息的通信。图5(b)和(c)分别描述了两个阶段的通信。在每一阶段，进程都会处理传递它在以前阶段所接收到的且需要再次传递的粒子。对第一阶段而言，进程只需通信自己区域内的粒子，而在此后所有的阶段，进程所传递的信息不仅包括了自己区域内的粒子，而且还包括从其他进程接受的粒子信息。

现参考图6、图7和图8，图6为根据本发明实施例的二维模拟系统二维分割边界区域第一次通信示意图；图7为根据本发明实施例的二维模拟系统二维分割边界区域第二次通信示意图；图8为根据本发明实施例的三维模拟系统三维分割边界区域的通信示意图。尽管在参考文献8(“Jiangxin Lu，Jiayuan Zhang，Xiaowei Wang，Limin Wang，Wei Ge.Particuology，7，317-323，2009”)中有类似的二维处理方式，但其通信方式在处理边界通信时通信的内容过多，而且计算一段时间后计算的效率会下降，需要定期对内存重构。在本发明的实施例中所提出的通信和处理方式则避免了这些情况，保证了计算的效率。

其中，如图6、7所示在二维模拟系统中，对于边界区的元胞网格，根据所相邻的子区域的个数可以区分为：只与一个子区域相邻的边界区元胞网格，也可以称为“棱”区域元胞网格；以及同时与两个子区域相邻的边界区元胞网格，也可以称为“顶点”区域元胞网格。而在如图8所示的三维模拟系统中将只与一个子区域相邻的边界区元胞网格称为“面心”区域元胞网格；将同时与三个子区域相邻的边界区元胞网格称为“棱”区域元胞网格；将同时与七个子区域相邻的边界区网格称为“顶点”区域元胞网格。处理边界区粒子碰撞时，首先处理相邻进程最少的边界区元胞网格，接着处理相邻进程次少的边界区元胞网格，最后才处理相邻进程最多的边界区元胞网格，依次顺序处理所有边界区元胞网格。在二维情形下，首先向高位进程传送高位“棱”区域中的粒子信息，在三维情形则首先传递“面心”区域网格中的粒子信息，处理完之后再传递“棱”区域中的粒子信息。其次，每个进程处理低位“棱”区域中粒子和从低位进程接收到的粒子之间的碰撞；碰撞完之后更新粒子的速度信息，并将从低位进程接收到的粒子更新后的信息反馈给低位进程；之后再处理“顶点”区域粒子的碰撞，并将更新后的粒子信息反馈给低位进程。这样对于二维系统二维分割情形边界区的粒子通信需要两次完整的通信处理，而对于三维系统三维分割情形则需要三次完整的通信处理，这种通信处理的增加是由于硬球模型对于碰撞处理的一致性所导致的。和参考文献8中不同的是，在通信处理的时候，为最大限度利用Shift通信模式和节省通信时间，那些对当次传递所不需要的粒子信息则不必参与传递，在传递顶点的网格中粒子的时候，顶点元胞网格的选取应该是分散的；并且在碰撞处理上，为了绝对确保碰撞的一致性，对邻近进程之间的粒子的碰撞处理只由一个进程处理，处理完之后再将信息反馈至相邻进程

针对在本实施例中的二维模拟系统，依次通过下列步骤来完成边界区“棱”区域元胞网格和“顶点”区域元胞网格中粒子的信息传递、碰撞处理以及信息反馈，具体步骤包括：

(1)利用计算节点计算进程进行边界区粒子的信息传递，实现低位进程的边界区高位“棱”区域元胞网格(图6中左斜线网格)中的粒子信息传递至邻近高位进程的虚拟元胞网格中(高位和低位进程参考图2)；

(2)计算子区域边界区中低位“棱”区域元胞网格中粒子和(1)通信后存储在虚拟区元胞网格(图6中右斜线网格)中粒子之间的碰撞，并更新碰撞粒子的状态；

(3)将已碰撞且更新状态后的虚拟区(图6中右斜线网格)的粒子信息，通过相邻进程间的通信反馈至邻近低位计算进程边界区对应的元胞网格之中，并从虚拟区中删除粒子的信息；

(4)将不同进程间的边界区高位的“顶点”区域元胞网格(图7中左斜线网格)中的粒子信息，以Shift通信方式传递至邻近进程虚拟元胞网格中；

(5)计算所述子区域边界区中低位“顶点”区域元胞网格(图7中右斜线网格)中，相互间尚未完成碰撞的粒子之间的碰撞，并更新碰撞粒子的状态；

(6)将所述已碰撞且更新状态后(图7中右斜线网格中)的粒子信息，通过不同计算进程间的Shift通信方式传递反馈至邻近计算进程边界区元胞网格中；

需要注意的是，如果模拟的系统是三维系统且并行算法为三维分割，则首先必须先处理“面心”区域(图8中所示)元胞网格中粒子的信息传递、碰撞处理以及信息反馈，其过程和(1)-(3)类似。

通过上述的操作，完成了子区域部分粒子时间驱动和事件驱动耦合计算和状态更新。

重复执行步骤1和步骤2，直至满足结束条件时将模拟后的结果从所述并行计算系统中输出。在一个实施例中，结束条件可以是指整个粒子系统达到稳定状态。在又一个实施例中，结束条件可以是指满足用户所要求的重复计算次数或计算时间。在其他实施例中，结束条件可以是指本领域普通技术人员所知的满足基于硬球模拟的其他结束条件，例如硬球之间的碰撞次数等。

一般而言，事驱区域应远大于时驱区域，在时驱区域，相当于空间换时间，而在事驱区域则相当于时间换空间，耦合算法同样也是对时空的耦合。

图9给出了基于硬球模型的大规模并行计算方法的一个实施例的流程图。最初，步骤a将所要模拟的粒子系统信息和模拟参数读入所述的并行计算系统中；模拟参数包括模拟体系的大小、所需进程数、进程维度分布以及一些用户所需参数等。步骤b，根据用户输入的模拟参数实现并行任务的划分，将所要模拟的系统划分为多个子区域；同时将所述子区域划分为多个元胞网格，所述子区域元胞网格分为边界区元胞网格和内部区元胞网格，同时在所述子区域外层建立子区域的虚拟区元胞网格；其中，内部区元胞网格对应事件驱动区域，边界区元胞网格和虚拟区元胞网格对应时间驱动区域。在又一个实施例中，也可以根据并行系统的进程总数实现并行任务的划分。步骤c，建立所要模拟的系统中粒子与元胞网格的映射关系，并完成事驱区域的预处理。步骤d，更新事件驱动区域内粒子的时间、速度和位置。步骤e，更新时间驱动区域粒子的时间和位置，并完成进程间粒子的迁移。步骤f，处理进程内时间驱动区域粒子的碰撞以及时间驱动区域粒子和事件驱动区域内粒子之间的碰撞。步骤g，处理邻近进程间粒子之间的碰撞并将更新后的粒子信息反馈给邻近进程。

重复执行步骤d、e、f、g，直至整个粒子系统达到稳定状态或满足用户所要求的重复计算次数或计算时间，将模拟后的结果从所述并行计算系统中输出。

以上就是时间驱动和事件驱动耦合实现和边界区粒子信息通信的实现步骤。采用上述方法可以得到良好的计算效果。

图10给出了单进程耦合算法和时驱算法以及事驱算法的效率比较，比较每步计算时间，以秒为单位。模拟体系的粒子数为1.6×10⁵，粒子的无量纲质量和半径均为1，填充率为0.0873，热速度为0.2，时驱算法和耦合算法的时间步长为1.0。图10示出了两个耦合事驱和时驱的算法：90％事件驱动区域和10％的时间驱动区域；50％事件驱动区域和50％的时间驱动区域。通过测算，该耦合算法的计算效率介于时间驱动和事件驱动算法之间，并随着时驱区域和事驱区域的变化而变化。可见，本发明实施例中所提供的耦合算法克服了硬球模型事件驱动算法难以并行且现有并行计算方法由于反复回调而导致的低效困难，而且其计算的效率要显著优于时间驱动算法。

图11是耦合并行算法和事驱单机算法对泊肃叶流动中流体的平均速度随时间的变化图。模拟中，粒子的无量纲质量和半径均为1，体系的宽度为1200，长度为1800，填充率为0.0873，外力场强为3.0×10^-8，初始温度为1.0，k_B为2.0×10^-4，耦合算法的时间步长为1.0。算例中有6×10⁴个粒子，耦合算法分为2×3个计算进程处理。从图中可以看出，对于泊肃叶流动中流体的平均速度随时间的演化，本发明实施例所提供的耦合事驱和时驱的并行方法与采用单机事驱算法模拟结果基本上相吻合。由此可见耦合事驱和时驱的并行方法的计算过程是准确和有效的，这同时也证实了本发明实施例提供的耦合事驱和时驱的并行方法的可靠性。

从上述试验可知，本发明实施例中所提供的基于硬球模型的并行计算方法不仅克服了硬球模型事件驱动算法难以并行且现有并行计算方法由于反复回调而导致的低效困难，同时也解决了时驱算法由于计算量大而造成的计算效率低下问题，在计算的效率和并行扩展性上都有了很大的提高。

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

Claims (11)

1.一种基于硬球模型的并行计算方法，所述方法包括以下步骤：

步骤1)将所要模拟的粒子系统信息读入并行计算系统中；

步骤2)将所要模拟的系统划分为多个子区域，同时将所述子区域划分为多个元胞网格并且在所述子区域外层建立子区域的虚拟区元胞网格，所述子区域元胞网格分为边界区元胞网格和内部区元胞网格；其中，所述内部区元胞网格对应事件驱动区域，所述边界区元胞网格和虚拟区元胞网格对应时间驱动区域；

步骤3)建立所要模拟的系统中粒子与元胞网格的映射关系；

步骤4)在给定的时间驱动时间步中，更新所述事件驱动区域和时间驱动区域内粒子的状态。

2.根据权利要求书1所述的并行计算方法，其特征在于，在所述步骤2)中，每个子区域分别对应所述并行计算系统中的一个进程。

3.根据权利要求书1所述的并行计算方法，其特征在于，在所述步骤2)中，所述虚拟区元胞网格和边界区元胞网格的层数均不少于2层，并且虚拟区元胞网格的层数不超过边界区元胞网格的层数。

4.根据权利要求书1所述的并行计算方法，其特征在于所述步骤3)中建立所要模拟系统中粒子与元胞网格的映射关系包括：根据粒子的位置将粒子以链表方式分布至步骤1)中所述子区域元胞网格中。

5.根据权利要求书1所述的并行计算方法，其特征在于，所述步骤4)中所述给定的时间驱动时间步Δt的大小应满足：

|v₁-v₂|Δt≤min(d/2，l₁₂-r₁-r₂)

其中，v₁、v₂是两个粒子的速度，l₁₂是粒子之间的距离，d是粒子的平均直径，r₁和r₂是粒子的半径。

6.根据权利要求书5所述的并行计算方法，其特征在于所述步骤4)中更新所述事件驱动区域内粒子的状态包括以下步骤：

步骤4-1)建立所述子区域内部区元胞网格的事件网格；

步骤4-2)将所述子区域中每个事件网格中的粒子的事件按发生的先后顺序存储在事件网格的一个链表中；

步骤4-3)找到所述子区域内部区元胞网格所对应的所有事件网格中最优先发生事件；

步骤4-4)如果所述最优先发生事件为穿越元胞网格事件则更新其时间和位置，并更新粒子所处元胞网格所在的事件网格的事件链表；同时分析粒子新位置所处区域，如果粒子处于时驱区域则更新粒子时间至时驱时间步并继续更新粒子的位置，如粒子的最新位置在事驱区域，则更新该最新位置所在元胞网格的事件网格的事件链表；

步骤4-5)如果所述最优先发生事件为碰撞事件则更新碰撞粒子对的时间、位置和速度；分析碰撞粒子对新位置所处元胞网格，并更新粒子链表发生更改的元胞网格所处事件网格中的事件链表；

步骤4-6)重复执行步骤4-3至4-5，直到处理完给定时间步内事驱区域所有事件。

7.根据权利要求书5所述的并行计算方法，其特征在于所述步骤4)中更新所述时间驱动区域内粒子的状态包括以下步骤：

步骤4-7)更新所述边界区元胞网格中的粒子的位置信息；

步骤4-8)、处理所述边界区元胞网格中粒子之间的碰撞，以及边界区元胞网格中粒子和内部区元胞网格中粒子的碰撞，并更新该粒子和被碰撞粒子的速度；

步骤4-9)、处理本子区域边界区元胞网格中粒子与相邻子区域边界区元胞网格中粒子之间的碰撞，并更新碰撞粒子对的速度。

8.根据权利要求书7所述的并行计算方法，其特征在于，所述步骤4-7)中更新粒子的位置信息包括以下步骤：

步骤4-7-1)、如果粒子的位置仍然处于原先的元胞网格，则粒子所属网格信息不变；

步骤4-7-2)、如果粒子的位置不在原先的元胞网格，但粒子的新位置仍处于所述边界区和虚拟区元胞网格中，则将粒子移至新位置所处的元胞网格中；如果粒子的位置不在原先的元胞网格，且粒子的新位置处于内部区元胞网格中，则将粒子移至新位置所处的元胞网格中，并更新新位置所处的元胞网格所在的事件网格的事件列表；

步骤4-7-3)、将虚拟区元胞网格中的粒子通信至相邻子区域，并从虚拟区元胞网格中删除粒子的信息。

9.根据权利要求书7所述的并行计算方法，其特征在于，所述步骤4-9)是依次按以下顺序处理处理本子区域边界区元胞网格中粒子与位于其它子区域边界区元胞网格中粒子之间的碰撞的：首先处理具有最少相邻子区域的边界区元胞网格，接着处理具有次少相邻子区域的边界区元胞网格，直到处理完具有最多相邻子区域的边界区元胞网格为止。

10.根据权利要求书9所述的并行计算方法，其特征在于，所述步骤4-9)处理本子区域边界区元胞网格中粒子与相邻子区域边界区元胞网格中粒子之间的碰撞包括以下步骤：

a.将本子区域的与高位子区域相邻的边界区元胞网格粒子的信息传递至所述高位子区域的虚拟区元胞网格中；

b.本子区域对应的进程计算与低位子区域相邻的边界区元胞网格中粒子和虚拟区元胞网格中粒子之间的碰撞，并更新碰撞粒子的状态；

c.将已碰撞且更新状态后的虚拟区的粒子信息反馈至对应的低位子区域边界区元胞网格之中，并从虚拟区中删除粒子的信息；

其中，所述高位子区域与低位子区域是按照步骤2)中所划分的子区域的位置而确定的。

11.根据权利要求1所述的并行计算方法，其中可以重复执行步骤4)直至系统达到稳定状态或者满足用户要求。

Images