CN110211235B - 基于并行rcb三维势函数离散元的放矿仿真方法 - Google Patents
基于并行rcb三维势函数离散元的放矿仿真方法 Download PDFInfo
- Publication number
- CN110211235B CN110211235B CN201910397522.2A CN201910397522A CN110211235B CN 110211235 B CN110211235 B CN 110211235B CN 201910397522 A CN201910397522 A CN 201910397522A CN 110211235 B CN110211235 B CN 110211235B
- Authority
- CN
- China
- Prior art keywords
- partition
- partitions
- old
- icpu
- list
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 161
- 238000004088 simulation Methods 0.000 title claims abstract description 24
- 230000008569 process Effects 0.000 claims abstract description 79
- 230000006854 communication Effects 0.000 claims abstract description 72
- 238000004891 communication Methods 0.000 claims abstract description 69
- 230000001133 acceleration Effects 0.000 claims abstract description 22
- 238000004364 calculation method Methods 0.000 claims abstract description 18
- 230000033001 locomotion Effects 0.000 claims abstract description 16
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 13
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 10
- 238000005457 optimization Methods 0.000 claims abstract description 6
- 238000005192 partition Methods 0.000 claims description 448
- 238000000638 solvent extraction Methods 0.000 claims description 16
- 238000003491 array Methods 0.000 claims description 7
- 125000004122 cyclic group Chemical group 0.000 claims description 7
- 238000013016 damping Methods 0.000 claims description 6
- 238000012856 packing Methods 0.000 claims description 5
- 238000001514 detection method Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 239000011435 rock Substances 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 abstract description 7
- 230000003993 interaction Effects 0.000 abstract description 6
- 238000006073 displacement reaction Methods 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 13
- 238000011439 discrete element method Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 2
- 230000005484 gravity Effects 0.000 description 2
- 108020001568 subdomains Proteins 0.000 description 2
- 101100450709 Mus musculus Hif3a gene Proteins 0.000 description 1
- 240000003864 Ulex europaeus Species 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2210/00—Indexing scheme for image generation or computer graphics
- G06T2210/32—Image data format
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了基于并行RCB三维势函数离散元的放矿仿真方法，通过采用RCB方法进行区域分解，以分配各进程的计算单元；采用MPI通信函数在各区域间进行通信；采用图论的方式进行区域间通信的免死锁优化；采用动态区域分解的方法保证各进程的负载平衡；计算过程中由牛顿第二定律求解颗粒的加速度与角加速度；通过velocity verlet算法计算块体的速度与位移。本发明解决了离散元并行计算中区域边界单元频繁运动引起进程间复杂交互的问题，解决了由于颗粒运动导致的负载不平衡问题，可用于在多核处理器或计算机集群上的多进程并行DEM仿真。
Description
技术领域
本发明涉及基于并行RCB三维势函数离散元的放矿仿真方法，可以实现超大规模三维非连续介质的变形与运动计算，属于高性能并行计算技术领域。
背景技术
离散单元法(Discrete Element Method，DEM)是一种拉格朗日方法，其根据牛顿第二运动定律计算单个粒子，从而在粒子水平上研究粒状流动特性，并精确地评估粒子运动。另一方面，当粒子数量很大时，计算成本变得过大，从而将DEM限制在小规模系统中，因此迫切需要开发新颖的建模和并行计算技术，以便在复杂和大规模的仿真过程中更好地进行DEM计算。
目前关于DEM的并行计算方法通常基于分布式内存或共享式内存环境开发，其中分布式内存环境多见于超级计算机和计算机集群中，相对于共享式内存环境而言，虽然其需要额外的通信开销，但其可扩展性强，可以处理规模更大的问题，且由于其每个进程都拥有自己的独立变量和内存，内存读写时不需要对数据进行频繁的加锁解锁。
目前大部分的分布式并行化策略通常试图将大型问题(计算域)划分为多个较小的子问题，各子问题由单独的子进程执行，各子进程之间则通过通信进行协调，离散元也是如此。其划分方法通常分为空间分解法、块体分解法和作用力分解法。其中，空间分解法把整个物理三维空间分解成若干个小的子空间并分配给每个进程，所以彼此的计算相对独立，克服了块体分解法中通信量非常大，和作用力分解法中重复计算的问题，因此得到的使用更加广泛。
在空间分解法中，Munjiza教授提出了一种基于二维的有限离散单元的并行算法，该算法通过RCB(Recursive Coordinate Bisection，递归坐标二分方法)进行区域分解，并在每个区域四周建立缓冲区，子域中由单元的坐标为单元划分类型，以此为依据在每个时步内进行相邻子域间单元的通信。但该方法的重叠区宽度很小，使得交互很复杂。事实上，重叠区的大小毕竟是降维的，占比很小，重叠区的存在使计算量的增加并不多，其次，该方法只解决了二维离散元的问题，对于三维空间的离散元问题其交互情况会变得更加复杂，而难以实现。
综合目前的研究，需要提出一种基于MPI并行计算的新的算法，去解决离散元并行计算中区域边界单元频繁运动引起进程间复杂交互的问题，以及由于块体运动导致的负载不平衡问题，可用于在多核处理器或计算机集群上的多进程并行DEM仿真。
发明内容
本发明所要解决的技术问题是：提供基于并行RCB三维势函数离散元的放矿仿真方法，针对目前离散元计算规模越来越大的情况，突破串行计算对于仿真规模的限制，并解决三维并行块体离散元中的进程间交互死锁的问题。
本发明为解决上述技术问题采用以下技术方案：
基于并行RCB三维势函数离散元的放矿仿真方法，包括如下步骤：
步骤1，为放矿建立块体离散单元模型，得到所有的块体单元，并对块体单元进行编号，各进程分别读取所有块体单元的位置坐标、几何信息、弹性模量、粘度、摩擦系数，全局仿真环境中的边界信息、重力加速度，以及边界的弹性模量、粘度、摩擦系数；
步骤2，确定分区的重叠区宽度，并采用RCB方法对步骤1的模型进行初次空间区域分解，得到分区，于此同时各进程对各块体单元进行分类标记，确定各块体单元位于分区的外重叠区、内重叠区或内区，并确定分区之间的关联列表，采用图论方法进行初始分区通信优化，确定初始分区通信的顺序和步骤；
步骤3，确定计算时间以及计算的时间步长，得到所有时间步；
步骤4，由上一个时间步计算所得的各块体单元的坐标、速度和加速度，预估当前时间步块体单元的速度，更新各块体单元的坐标；
步骤5，判断是否需要重新进行RCB分区，如果需要，则重新RCB分区，记为新分区，并将步骤2的分区作为老分区，将老分区的单元列表和形心、速度、角速度物理量传递给与之关联的新分区；否则进入步骤6；若重新RCB分区，则步骤6-9所述分区即为新分区，否则即为老分区；
步骤6，根据步骤4更新的坐标，各进程对穿越分区的块体单元进行信息交换；
步骤7，各进程对分区内的块体单元进行接触检测，得到接触链表，与此同时计算当前时间步对应分区内且包含外重叠区每个块体单元所受的接触力和力矩；
步骤8，将各分区内重叠区中块体单元接触力和对应相邻分区的外重叠区中块体单元接触力进行交换；
步骤9，由刚体运动方程，根据接触力计算出每个块体单元当前时间步的加速度，再由verlet算法计算出每个块体单元下一时间步的速度；
步骤10，重复步骤4至步骤9计算下一时间步，直至计算完所有时间步，得到放矿后岩石散落的形态。
作为本发明的一种优选方案，所述步骤2的具体过程如下：
2.1计算各块体单元的形心坐标，求块体单元各角点与形心之间的距离，取最大距离作为块体单元的半径，取半径的2.01倍作为块体单元的直径，将所有块体单元对应的直径中最大直径的1.05倍作为内、外重叠区的宽度；
2.2根据块体单元的编号，将块体单元平均分配到各进程中；
2.3创建保存初始分区信息的自定义类型cpuRCB_type，类型中包含四个变量，分别为：分区单元数量ne、单元索引列表liste、分区左下角坐标Domleftlow、分区尺寸Domsize；
2.3.1cpuRCB_type类型数组长度的确定
指定每个分区方向的分区层数，由RCB二分的性质得到：
(1)若第一个分区方向的分区层数是N1，则该方向最后一层分区条带数
总分区数
(2)若第二个分区方向的分区层数是N2，则该方向最后一层分区条带数
总分区数
(3)若第二个分区方向的分区层数是N3，则该方向最后一层分区条带数
总分区数
令N＝N1+N2+N3，则最后一层总分区数ncpu＝M1×M2×M3＝2N-3；
所有层的总分区数ncpuRCB＝L1×L2×L3，以ncpuRCB为大小创建cpuRCB_type类型数组；
2.3.2变量ne、liste、Domleftlow和Domsize的确定
(1)设定分区单元数量允许容差ratioRCB和二分法最大迭代次数niterRCB，指定二分方向为idimn，idimn＝1,2,3分别表示二分方向为X,Y,Z；
(2)以指定idimn方向的空间区域两侧坐标为边界，以中间坐标为界限将空间区域划分为两个区域，以二分方向idimn＝1为例：初始空间区域左下角顶点坐标为(X0,Y0,Z0)，空间区域在三个方向上的长度分别为[Lx,Ly,Lz]，则该方向以X1＝X0+Lx/2作为分割线，分割为两个区域：
区域1：左下角坐标为(X0,Y0,Z0)，区域尺寸为[X1-X0,Ly,Lz]；
区域2：左下角坐标为(X1,Y0,Z0)，区域尺寸为[X0+Lx-X1,Ly,Lz]；
当块体单元的形心(Cx,Cy,Cz)都在区域1内时，即：
X0<Cx<X1-X0
Y0<Cy<Ly
Z0<Cz<Lz
若上述三个条件都满足，则判断该块体单位在区域1内，同理判断单元是否在区域2内，由此统计两个区域内的分区单元数量ne和单元索引列表liste；
(3)各进程归约统计得到区域1和2块体单元的实际总数量，若区域1块体单元总数量大于区域2块体单元总数量，则令新的分割线为X2＝(X0+X1)/2，反之，则令新的分割线为X2＝(X0+X1+Lx)/2；
(4)重复步骤(1)-(2)，直到区域1和2块体单元总数量的差值与总块体单元数量的比值小于允许容差ratioRCB时，或者循环次数大于最大迭代次数niterRCB时，退出循环，此时得到的分割线坐标Xn即为二分法要得到的值，且二分法后两个区域的空间位置：
区域1：左下角坐标为(X0,Y0,Z0)，区域尺寸为[Xn-X0,Ly,Lz]；
区域2：左下角坐标为(Xn,Y0,Z0)，区域尺寸为[X0+Lx-Xn,Ly,Lz]；
(5)将区域1的左下角坐标、尺寸分别存入cpuRCB_type类型数组的第1个元素中，区域2的左下角坐标、尺寸分别存入cpuRCB_type类型数组的第2个元素中；
(6)将上述(2)-(4)的操作称为idimn方向的第1层分区，则第2层分区为：将第1层分区得到的两个区域分别作为(2)中的初始空间区域，重复(2)-(3)，得到的4个子分区，存入cpuRCB_type类型数组的第3至第6个元素中；以此类推，直至idimn方向执行完第Nidimn层分区；
(7)按照上述(6)中的方法进行其他两个方向的分区，并填充到cpuRCB_type类型数组中；
2.4建立分区关联列表
2.4.1获取分区三维编号
分区三维编号由其所在空间分区拓扑中的位置确定，例如，分区编号(I,J,K)代表该分区在所有分区中的排序为：第I行(X方向)，第J列(Y方向)，第K层(Z方向)；
2.4.2确定进程编号和分区编号的关系
通过MPI_CART_CREATE函数创建笛卡尔拓扑通信域Cart_Comm，笛卡尔拓扑的各方向进程数量与分区在各方向的数量一致，由此，在通信域Cart_Comm中，由MPI_CART_RANK函数获取分区三维编号(I,J,K)所对应的进程编号，通过MPI_CART_COORDS函数由进程编号获取分区三维编号(I,J,K)；
2.4.3获取各分区的内外重叠区
将2.3中划分好的各分区的边界E向内收缩一个重叠区宽度，得到边界E内，将边界E向外扩展一个重叠区宽度得到边界E外，则边界E外到边界E之间的区域称为该分区的外重叠区，边界E到E内之间的区域称为该分区的内重叠区，边界E内以内的区域称为该分区的内区；
根据各分区中块体单元所在位置将块体单元标记为如下三种情况：
(1)若块体单元位于内区，标记为1；
(2)若块体单元位于内重叠区，标记为2；
(3)若块体单元位于外重叠区，标记为3；
2.4.4分区关联判断具体步骤：
(1)创建cpuDB_type类型的数组，数量为最后一层总分区数ncpu，cpuDB_type类型包含分区左下角坐标、分区尺寸、关联分区的个数ndlinkDB、关联分区列表dlinkDB；将2.3中cpuRCB_type类型数组第2N-1到2N-1个元素中分区的几何信息分别复制到cpuDB_type类型数组的第1到2N-1个元素中；
(2)确定分区三维编号到cpuDB_type类型数组的映射关系
分区的三维编号(I,J,K)与cpuDB_type类型数组中的编号icpu关系如下：
第I个面、第J个条带的第K个分区编号icpu是：
icpu＝(I-1)×M2×M3+(J-1)×M3+K
分区数组编号为icpu的分区，其三维编号(I,J,K)是：
I＝(icpu-1)/(M2×M3)+1
J＝(icpu-(I-1)×M2×M3-1)/M3+1
K＝icpu-(I-1)×M2×M3-(J-1)×M3
(3)筛选出有可能关联的分区
对三维编号为(I,J,K)的分区来说，有可能关联的分区为如下两类：
①一定会通讯的是：(I,J,K+1)
②可能会通讯的是：(I,J+1,1:M3)，(I+1,1:M2,1:M3)
其中，1:M3表示从1到M3的数组，1:M2表示从1到M2的数组；
(4)对可能关联的分区进行关联判断
对于cpuDB_type类型数组中编号为icpu的分区而言，要对所有可能关联的分区进行判别，具体判别方法如下：
①将可能通讯的分区三维编号(Ij,Jj,Kj)转化为cpuDB_type数组中编号jcpu；
②获取cpuDB_type数组中编号为icpu和jcpu的元素中的分区几何信息，即左下角坐标与尺寸，将左下角坐标的三个值分别加上三个方向的尺寸，求得分区三个方向的范围；对于分区序列(I,J+1,1:M3)中的分区，判断第三个方向是否在(I,J,K)分区第三个方向范围内，若在则分区关联，否则不关联；对于分区序列(I+1,1:M2,1:M3)中的分区，判断第二和第三个方向是否在(I,J,K)分区第二和第三个方向范围内，只要有一个方向不在范围内，则分区不关联；
2.5优化通讯批次
采用图论的方法优化通讯的顺序批次，保证每个批次中没有重复的通讯进程，具体步骤如下：
(1)标记cpuDB_type数组每个元素中所有关联分区为未使用状态，再将所有分区标记为未循环状态；
(2)从cpuDB_type数组中的第一个元素开始索引，如果该分区为未循环状态，则寻找该元素关联分区列表中第一个处于未使用且未循环状态的分区序号，如果有，则将其序号jcpu与该元素序号icpu记为一个通信对，最后将jcpu标记为循环且已使用状态；重复该操作，直到循环完所有元素，则完成了一批通讯批次；
(3)重复第(2)步，直到在某个批次中找不到通信对为止；
(4)按照求得的通讯批次的先后顺序，将cpuDB_type数组每个元素对应的关联分区编号重新填写；
(5)每一轮通信时，编号小的进程先接受后发送，编号大的进程先发送后接受。
作为本发明的一种优选方案，步骤3所述时间步长小于等于允许的最大时间步长，且允许的最大时间步长为：
其中，Δt是允许的最大时间步长，ξ是阻尼比，
c是阻尼系数，m是块体单元的质量，k是刚度系数。
作为本发明的一种优选方案，所述步骤4的具体过程如下：
当前时间步块体单元的速度，预估公式为：
更新坐标由下式进行：
其中，
为当前时间步预估速度，
为上一时间步速度，Δt为时间步长，
为上一时间步加速度，un+1为当前时间步位移，un为上一时间步位移。
作为本发明的一种优选方案，所述步骤5的具体过程如下：
5.1判断是否需要进行重新RCB分区的标准如下：
(1)根据总的块体单元数/总分区数ncpu，得到平均单元数；
(2)对当前分区，若块体单元数小于平均单元数的某个比例，则认为该分区是坏的，赋值标记为1，否则为0；
(3)归约统计坏分区的数量；
(4)如果坏分区的数量大于总分区数ncpu的某个比例，则需要重新进行RCB分区；
5.2重新RCB方法：
5.2.1重复步骤2.3-2.4完成对空间区域的重新划分，得到新分区；
5.2.2新老分区关联
(1)判断是否关联
将每一个老分区与每一个新分区在三个方向比较其范围，若老分区在某个方向坐标的最小值大于新分区对应方向坐标的最大值，或在某个方向坐标的最大值小于新分区对应方向坐标的最小值，则说明该老分区与新分区无关联，否则有关联；
(2)将每个老分区关联的新分区编号记录到新老分区关联列表中；
5.2.3优化新老分区间的通讯批次
(1)分离新老分区关联列表：
对icpu号老分区的关联列表进行索引，若关联列表中有jcpu号新分区，再索引jcpu号新分区的关联列表，若存在icpu号老分区，则在jcpu号老分区的关联列表中剔除icpu号，并在新建的关联列表中jcpu号老分区关联的队列中添加-icpu号；
(2)将步骤2得到的两个关联列表调整为只包含比自身编号大的分区的状态：
icpu号老分区将其关联列表中的所有的分区标记为关联状态，再索引比icpu号大的老分区关联列表中是否有icpu号新分区，若有则将icpu号老分区标记关联状态，并将比icpu号大的老分区关联列表中的icpu号剔除，最后将所有老分区中处于关联状态的进程从小到大写入icpu号老分区关联列表中；由此方法按编号从小到大依次重填关联列表，即可得到关联列表只包含比自身编号大的分区的列表；
(3)重复步骤2.5，将上一步中所得的两个新关联列表分别优化，得到两个通讯批次；
(4)合并两个新关联列表，得到合并后的关联列表；
(5)按照5.2.2得到的新老分区关联列表标记合并后的关联列表中的发送和接受状态；
5.2.4新老分区单元信息互换
(1)创建dlinkDB_type类型数组dlinkDB，该数组长度为本老分区关联分区的数量，类型包含5个变量：关联的分区编号jcpu、需要传送给jcpu的单元索引列表listout、传送的单元个数neout、穿越到jcpu的单元列表listpass、穿越到jcpu的单元数量nepass；
(2)统计各分区内的块体单元数量及单元索引列表，分别存储到dlinkDB数组中的变量neout和listout中；
(3)按照5.2.3中合并后的通讯批次，按照dlinkDB中neout发送单元数量；
(4)按照上述(3)中的通讯批次，交换listout数组；
(5)按照索引列表listout将需要发送的单元的形心、速度、角速度打包到一个数组A中；
(6)将数组A发送到分区jcpu中，分区jcpu接受后，将接受到的数组解包，并填充到对应单元的形心、速度、角速度列表中。
作为本发明的一种优选方案，所述步骤9的具体过程如下：
计算加速度的公式如下：
计算下一时间步速度的公式如下：
其中，
为当前时间步加速度，Fn+1为当前时间步接触力，m为块体单元质量，
为下一时间步速度，
为当前时间步预估速度，Δt为时间步长。
本发明采用以上技术方案与现有技术相比，具有以下技术效果：
本发明解决了离散元并行计算中区域边界单元频繁运动引起进程间复杂交互的问题，解决了由于颗粒运动导致的负载不平衡问题，可用于在多核处理器或计算机集群上的多进程并行DEM仿真。
附图说明
图1是本发明基于并行RCB三维势函数离散元的放矿仿真方法的流程图。
图2是RCB分区。
图3是块体单元跨越分区示意图。
图4是新老分区关联示意图。
图5是基于图论的分区通讯优化示意图。
图6是实施例放矿对应不同进程下的计算时间示意图。
图7是实施例放矿的块体离散单元模型。
图8(a)是放矿前的4进程分区的划分。
图8(b)是4进程下放矿时块体的运动过程。
图8(c)是4进程下放矿时的动态分区过程。
图9(a)是放矿前的8进程分区的划分。
图9(b)是8进程下放矿时块体的运动过程。
图9(c)是8进程下放矿时的动态分区过程。
具体实施方式
下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。
如图1所示，本发明提供了基于并行RCB三维势函数离散元的放矿仿真方法，包括以下步骤：
步骤1、各进程分别读取所有块体单元位置坐标、几何信息，以及全局仿真环境中的几何边界信息、重力加速度、块体及边界的弹性模量、粘度、摩擦系数等物理变量。具体为：
1.1各进程按单元序号分别读取相应序号区间内的块体单元角点坐标、几何信息、材料信息，以及全局仿真环境中的重力加速度等物理变量。
1.2对每个块体单元求解其势函数值、势多面体角点坐标、旋转变换矩阵、转动惯量、形心、最小外接球。
步骤2、确定重叠区宽度，并采用继承性RCB(Recursive Coordinate Bisection)方法进行初次空间区域分解，于此同时在各进程对各单元进行分类标记，并确定初始子分区关联列表，采用图论的方法进行子分区通信优化，确定分区通信的顺序和步骤。具体为：
2.1计算出各块体的形心坐标，将块体各角点与形心求距离，取最大距离作为块体的半径，取半径的2.01倍作为块体直径，将内外重叠区宽度均取为最大块体直径的1.05倍；
2.2按单元号顺序，将单元平均分配到各进程；
2.3创建保存RCB分区信息的自定义类型cpuRCB_type，类型中包含四个变量，分别为：该分区的单元数量ne、单元索引列表liste、分区左下角坐标Domleftlow、分区尺寸Domsize；
进一步的，针对“2.3”中cpuRCB_type类型数组的分配和求解，具体步骤如下：
2.3.1cpuRCB_type类型数组长度的确定
首先指定每个分区方向的分区层数，由RCB二分的性质可知每个方向第a层的分区数量为2a-1，若某个方向总分层数量为A，由等比数列求和公式可知每一层分区数之和ncpuRCB为：
具体的：
(1)若第一个分区方向的分区层数是N1，则该方向最后一层分区条带数
总分区数
(2)若第二个分区方向的分区层数是N2，则该方向最后一层分区条带数
总分区数
(3)若第二个分区方向的分区层数是N3，则该方向最后一层分区条带数
总分区数
令N＝N1+N2+N3，则最后一层总分区数
所有层的总分区数ncpuRCB＝L1×L2×L3，再以ncpuRCB为大小创建cpuRCB_type类型数组。
2.3.2变量Domleftlow、Domsize、ne和liste的确定
(1)设定分区单元数量允许容差ratioRCB和二分法迭代的最大层数niterRCB，指定二分方向为idimn(idimn＝1,2,3分别表示二分方向为X,Y,Z)。例如，指定二分方向的顺序为3、1、2，每个方向划分2层，则最终分区如图2所示。
(2)以指定idimn方向的分区两侧坐标为边界，以该分区的中间坐标为界限将分区划分为两块，例如：初始分区左下角顶点坐标为(X0,Y0,Z0)，区域三个长度分别为[Lx,Ly,Lz]，二分方向idimn＝1，则该方向则以X1＝X0+Lx/2作为分割线，可以分割为两个区域：
区域1：左下角坐标为(X0,Y0,Z0)，区域尺寸为[X1-X0,Ly,Lz]；
区域2：左下角坐标为(X1,Y0,Z0)，区域尺寸为[X0+Lx-X1,Ly,Lz]；
当块体单元的形心(Cx,Cy,Cz)都在区域1内时，即：
X0<Cx<X1-X0
Y0<Cy<Ly
Z0<Cz<Lz
三个条件都满足，同理判断单元是否在区域2内，由此统计两个子分区内的单元数量ne和单元列表liste。
(3)各进程归约统计得到两侧单元的实际总数量，若区域1单元数量大于区域2，则令新的分割线为X2＝(X0+X1)/2，反之，则令新的分割线为X2＝(X0+X1+Lx)/2；
(4)重复步骤(1)、(2)，直到最终两侧单元数量的差值与总单元数量的比值小于指定允许容差ratioRCB时，或者循环次数大于指定最大迭代层数niterRCB时，退出循环，此时得到的分割线坐标Xn即为二分法要得到的值，此时，可得到二分法后两个子分区的空间位置：
区域1：左下角坐标为(X0,Y0,Z0)，区域尺寸为[Xn-X0,Ly,Lz]；
区域2：左下角坐标为(Xn,Y0,Z0)，区域尺寸为[X0+Lx-Xn,Ly,Lz]；
(5)将区域1的左下角坐标、尺寸、分别存入cpuRCB_type类型数组的第1个元素中，区域2的左下角坐标与尺寸分别存入cpuRCB_type类型数组的第2个元素中。
(6)将上述“(2)”-“(4)”的操作称为idimn方向的第1层分区，则第2层分区可以叙述为：将第1层分区得到的两个子区分别作为“(2)”中的“初始分区”，重复“(2)”-“(3)”，得到的4个子分区，存入cpuRCB_type类型数组的第3到第6个元素中。
以此类推，直至idimn方向执行完第Nidimn层分区。
(7)按照上述“(6)”中的方法进行其他两个方向的分区，并填充到cpuRCB_type类型数组中。
2.4建立分区关联列表
2.4.1获取分区三维编号
分区三维编号由其所在空间分区拓扑中的位置确定，例如，分区编号(I,J,K)代表该分区在所有分区中的排序为：第I行(X方向)，第J列(Y方向)，第K层(Z方向)。
2.4.2确定进程编号和分区编号的关系
通过MPI_CART_CREATE函数创建笛卡尔拓扑通信域Cart_Comm，笛卡尔拓扑的各方向进程数量与分区在各方向的数量一致，由此，在通信域Cart_Comm中，可由MPI_CART_RANK函数获取分区三维编号(I,J,K)所对应的一维进程编号，通过MPI_CART_COORDS函数得到进程对应分区的三维编号(I,J,K)。
2.4.3获取子分区的内外重叠区
将“2.3”中划分好的各子分区的边界E向内缩一个重叠区宽度，得到边界E内，将边界E向外扩展一个重叠区宽度得到边界E外，则边界E外到边界E之间的区域称为该子分区的“外重叠区”，边界E到E内之间的区域称为“内重叠区”，边界E内所包围的区域称为“内区”。
根据各子分区中单元所在位置将单元标记为如下三种情况：
(1)内区，标记为1；
(2)内重叠区，标记为2；
(3)外重叠区，标记为3；
2.4.4分区关联判断具体步骤：
(1)创建cpuDB_type类型的数组，数量为分区总数量ncpu，cpuDB_type类型包含分区的几何信息(左下角坐标、三个方向尺寸)、关联分区的个数ndlinkDB、关联分区列表dlinkDB(:)。将步骤“2.3”中cpuRCB_type类型数组第2N-1到2N-1个元素(N为分区总层数)中分区的几何信息分别复制到cpuDB_type类型的数组中的第1到2N-1个元素中。
(2)确定分区三维编号到cpuDB_type类型数组的映射关系
分区的空间编号(I,J,K)与cpuDB_type类型数组中的编号icpu关系如下：
第I个面、第J个条带的第K个分区编号icpu是：
icpu＝(I-1)×M2×M3+(J-1)×M3+K
分区数组编号为icpu的分区，其空间编号(I,J,K)是：
I＝(icpu-1)/(M2×M3)+1
J＝(icpu-(I-1)×M2×M3-1)/M3+1
K＝icpu-(I-1)×M2×M3-(J-1)×M3
(3)筛选出有可能关联的分区
对空间编号为(I,J,K)的分区来说，有可能的关联是如下两类(遇到边界则为空)：
①一定会通讯的是：(I,J,K+1)
②可能会通讯的是：(I,J+1,1:M3)，(I+1,1:M2,1:M3)
其中，M2、M3为步骤“(2)”中的第二和第三个方向最后一层分区数，1:M3表示从1到M3的数组，1:M2表示从1到M2的数组；
(4)对可能关联的分区进行关联判断
对于cpuDB_type类型数组中编号为icpu的分区而言，要对所有可能关联的分区进行判别，具体判别方法如下：
①将可能通讯的分区空间编号(Ij,Jj,Kj)转化为cpuDB_type数组中编号jcpu；
②获取cpuDB_type数组中编号为icpu和jcpu的元素中的分区几何信息，即左下角坐标与尺寸，将左下角坐标的三个值分别加上三个方向的尺寸，可以求得分区三个方向的范围。对于分区序列(I,J+1,1:M3)中的分区，只需要判断第三个方向是否在(I,J,K)分区第三个方向范围内即可；对于分区序列(I+1,1:M2,1:M3)中的分区，则需要判断第二和第三个方向是否在(I,J,K)分区第二和第三个方向范围内，只要有一个方向不在范围内，则分区不关联。
2.5优化通讯批次
现将采用图论的方法优化通讯的顺序批次，以保证每个批次中没有重复的通讯进程，其具体步骤如下：
(1)标记cpuDB_type数组每个元素中所有关联分区为未使用状态，再将所有分区标记为未循环状态。
(2)从cpuDB_type数组中的第一个元素开始索引，如果该分区为未循环状态，则寻找该元素关联分区列表中第一个处于未使用且未循环状态的分区序号，如果有，则将其序号jcpu与该元素序号icpu记为一个通信对，最后将jcpu标记为循环且已使用状态。重复该操作，直到循环完所有元素，则完成了一批通讯批次。
(3)重复第(2)步，直到在某个批次中找不到通信对为止；
(4)按照求得的通讯批次的先后顺序，将cpuDB_type数组每个元素对应的关联分区编号重新填写，例如，现有6个进程，对应于6个分区，其空间关系如图5所示，其分为6个通讯批次如下所示：
第一批：(1，2)；(3，5)
第二批：(1，4)；(2，3)；(5，6)
第三批：(1，5)；(2，4)；(3，6)
第四批：(2，5)
第五批：(2，6)；(4，5)
则按照顺序，从第一批到第五批
与1号分区关联的为2、4、5、0、0
与2号分区关联的为1、3、4、5、6
与3号分区关联的为5、2、6、0、0
与4号分区关联的为0、1、2、0、5
与5号分区关联的为3、6、1、2、4
与6号分区关联的为0、5、3、0、2
其中，编号为0则表示该轮不通信。
(5)每一轮通信时，编号小的进程先接受后发送，编号大的进程先发送后接受，例如，在上例的第一轮通信中，1号与2号要互相通信，则1号先接受后发送，2号先发送后接受。
步骤3、确定计算时间步长。
允许的最大时间步长为：
其中，ξ是系统的阻尼比，
c是阻尼系数，m是块体单元质量，k是刚度系数。
步骤4、由上一个时步计算所得的各单元的坐标、速度和加速度，预估当前时步速度，更新各单元坐标。
由下式预估当前步速度：
式中，
为当前步预估速度，
为上一步速度，Δt为时间步长，
为上一步加速度。
更新坐标由下式进行：
式中，un+1为当前步位移，un为上一步位移。
步骤5、判断是否需要重新RCB分区，如果需要，则重新RCB分区，并交换新老分区之间单元的列表和物理量。具体为：
5.1重新RCB判断的标准如下：
(1)总的单元数/总分区数，得到的平均单元数；
(2)对当前分区，若单元数小于平均单元数的某个比例(0-0.1)，则认为这个分区是坏的，赋值标记为1，否则为0；
(3)归约统计坏分区的数量；
(4)如果坏分区的数量大于总分区数ncpu的某个比例(0.2-0.3)，则需要重新RCB。
5.2重新RCB方法：
5.2.1重复步骤“2.3”～“2.4”完成对区域的重新划分；
5.2.2新老分区关联
(1)判断是否关联
将每一个老分区与每一个新分区在三个方向比较其范围，若老分区在某个方向坐标的最小值大于新分区对应方向坐标的最大值，或在某个方向坐标的最大值小于新分区对应方向坐标的最小值，则说明该老分区与新分区无关联，否则说明有关联。例如：如图4所示，老分区中某个分区其左下角坐标为(1.0，1.0，1.0)，三个方向尺寸为(2.0，3.0，4.0)，则该分区在三个方向的范围分别为1.0～3.0、1.0～4.0、1.0～5.0，同理，若有某个新分区其三个方向范围为1.0～4.0、1.5～3.5、5.2～6.2，分别比较其三个方向范围可知，老分区第三个方向的最大值为5.0，而新分区第三个方向的最小值为5.2，大于5.0，因此该新老分区不关联。
(2)将每个老分区关联的新分区编号记录到数组D中,为方便下面说明，记Di[j]表示i号老分区在关联列表中的第j个关联的新分区的分区号。
5.2.3优化新老分区间的通讯批次
(1)分离新老分区关联列表：
对icpu号分区的关联列表进行索引，若关联列表中有jcpu号分区，再索引jcpu号分区的关联列表，若存在icpu号分区，则在jcpu号分区的关联列表中剔除icpu，并在新的关联列表中添加icpu。
(2)将关联列表调整为“关联列表只包含比自身编号大的分区”的状态：
icpu号老分区将其关联列表中的所有的分区进程标记为关联状态，再索引比icpu号大的老分区关联列表中是否有icpu号新分区，若有则将icpu号进程标记关联状态，并将关联列表中的icpu号剔除，最后将所有进程中处于关联状态的进程从小到大写入icpu号老分区关联列表中。由此方法按编号从小到大依次重填关联列表，即可得到“关联列表只包含比自身编号大的分区”的列表。
(3)重复步骤“2.5”，将上一步中所得的两个新关联列表分别优化，得到两个通讯轮次。
(4)合并两个关联列表。
(5)标记分区发送和接受，发送为正，接受为负，方法如下：
对于合并后的关联列表C，Ci[j]表示i号老分区在合并后的关联列表中的第j个关联的新分区的分区号：
一、则当j＝1时，i号老分区的第1个关联的新分区为Ci[1]，此时，遍历步骤“5.2.2”中第“(2)”步中的关联列表Di(即i号老分区在初始关联列表中的所有关联的新分区)，若其中存在某个新分区的编号与Ci[1]相同，则在关联列表C中，将Ci[1]标记为正，不存在则标记为负；
二、令j依次等于1、2…，直到关联列表的最后一项，并重复上述步骤“一”，即可完成对所有关联列表的标记。
5.2.4新老分区单元信息互换
(1)创建dlinkDB_type类型数组dlinkDB，该数组长度为本进程关联分区的数量。类型包含5个变量：关联的分区编号jcpu、需要传送给jcpu的单元索引列表listout、传送的单元个数neout、穿越到jcpu的单元列表listpass，穿越到jcpu的单元数量nepass(后两个变量在本步骤中不需要管理)；
(2)统计各分区内的单元数量及单元列表，分别存储到dlinkDB数组中的变量neout和listout中；
(3)按照“5.2.3”中合并后的通讯轮次，按照dlinkDB中neout发送单元数量；
(4)按照上述“(3)”中的通讯批次，交换listout数组；
(5)按照索引列表listout将需要发送的单元的形心、速度、角速度打包到一个数组A中；
(6)将数组A发送到进程jcpu中，对方进程接受后，将接受到的数组解包，并填充到对应单元的形心、速度、角速度列表中。
进一步的，上述“按照两批通讯批次”的具体方法如下：
为说明清楚，现举一例，以下为新老分区的关联列表(不包含与自身关联的新分区)：
与1号老分区关联的新分区为3、5、6
与2号老分区关联的新分区为1、3、4
与3号老分区关联的新分区为1、5、6
与4号老分区关联的新分区为2、6
与5号老分区关联的新分区为4、3
与6号老分区关联的新分区为5、1
分离关联列表如下：
调整为“关联列表只包含比自身编号大的分区”的状态如下：
分别优化后通讯轮次如下：
标记发送接受后，并合并的通讯轮次列表如下：
在上述列表中，负值表示接受，正值表示发送，每一轮按照一列进行通信，如第1轮通信为：1接受2；2发送给1；3发送给5；4发送给6；5发送至3；6接受4。
第二轮为：1发送给3；2发送给4；3接受1；4接受2；5接受6；6发送给5。
以此类推进行共7个批次的通讯。
步骤6、各进程对穿越分区的单元进行信息的交换。具体为：
6.1统计内重叠区内的单元数量neInBuf和单元列表listInBuf；
6.2将步骤“5.2.4”中dlinkDB数组各元素内的listout数组和listpass释放，并将neout和nepass置零；
6.3遍历单元列表listInBuf中的单元：若单元ie的当前步形心位于本进程的关联分区jcpu的外重叠区中，则将单元ie的索引添加到dlinkDB数组中分区编号为jcpu的元素中的listout数组中，并将neout数值加1；
6.4遍历dlinkDB数组每个元素：对于第id个元素，其分区编号为jcpu，遍历该元素中的listout数组，若单元ie形心上一步位于某个关联分区jcpu的外重叠区以外，则将该单元添加到listpass中，并将nepass数值加1。例如，在图3中，单元e当前步在B的外重叠区中，而上一步在B的外重叠区之外，且在A区之中，则分区A需要将单元e的索引添加到listpass中。
6.5按照“5.2.3”中优化好的通讯轮次，按照dlinkDB中nepass发送单元数量；
6.6按照通讯轮次，交换listpass数组；
6.7按照索引列表listpass将需要发送的单元的形心、速度、角速度打包到一个数组A中；
6.8将数组A发送到进程jcpu中，对方进程接受后，将接受到的数组解包，并填充到对应单元的形心、速度、角速度列表中。
步骤7、各进程对子域内的单元进行接触检测，得到接触链表，与此同时计算当前时间步对应子域内(包含外重叠区)每个块体单元所受的接触力和力矩。
步骤8、将各分区内重叠区和对应相邻分区的外重叠区单元接触力进行交换，使得各单元接触力信息准确。具体为：
8.1遍历dlinkDB数组：对于元素id，按照该元素索引列表listout将需要发送的单元的接触力和弯矩打包到一个数组F中；
8.2按照通讯轮次发送和接受数组F；
8.3将接受到的数组F解包，并赋值到对应各单元的接触力和力矩上。
步骤9、由刚体运动方程，根据接触力计算出每个块体单元当前时间步的加速度，再由verlet算法计算出每个块体单元下一时间步的速度。具体为：
计算加速度的公式如下：
式中，
为当前步加速度，Fn+1为当前步接触力，m为单元质量。
计算下一时间步速度的公式如下：
式中，
为下一步速度，
为当前步预估速度，Δt为时间步长。
步骤10、重复步骤4至步骤9计算下一时间步，直至计算完所有时间步。
实施例：
某矿区放矿，为研究崩落矿岩的流动特性与放矿规律，从而指定合理的放矿控制条件。采用本发明提供的方法，为放矿建立块体离散单元模型，如图7所示，共划分离散单元4415个。
基于本发明所提供的基于MPI的分布式三维块体离散单元并行算法，分别使用串行方法、4进程、8进程模拟放矿过程，各进程下计算所用时间如图6所示，通过比较计算时间可以发现，采用该方法可以对现有的离散元计算进行很好的加速效果。
图8(a)为放矿前的4进程分区的划分，图8(b)、图8(c)为4进程下放矿时块体的运动过程和动态分区过程。
图9(a)为放矿前的8进程分区的划分，图9(b)、图9(c)为8进程下下放矿时块体的运动过程和动态分区过程。
综上所述，采用基于MPI的分布式三维块体离散单元并行算法，能够对三维块体离散元的模拟提供良好的加速作用。
以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。
Claims (6)
1.基于并行RCB三维势函数离散元的放矿仿真方法，其特征在于，包括如下步骤：
步骤1，为放矿建立块体离散单元模型，得到所有的块体单元，并对块体单元进行编号，各进程分别读取所有块体单元的位置坐标、几何信息、弹性模量、粘度、摩擦系数，全局仿真环境中的边界信息、重力加速度，以及边界的弹性模量、粘度、摩擦系数；
步骤2，确定分区的重叠区宽度，并采用RCB方法对步骤1的模型进行初次空间区域分解，得到分区，于此同时各进程对各块体单元进行分类标记，确定各块体单元位于分区的外重叠区、内重叠区或内区，并确定分区之间的关联列表，采用图论方法进行初始分区通信优化，确定初始分区通信的顺序和步骤；
步骤3，确定计算时间以及计算的时间步长，得到所有时间步；
步骤4，由上一个时间步计算所得的各块体单元的坐标、速度和加速度，预估当前时间步块体单元的速度，更新各块体单元的坐标；
步骤5，判断是否需要重新进行RCB分区，如果需要，则重新RCB分区，记为新分区，并将步骤2的分区作为老分区，将老分区的单元列表和形心、速度、角速度物理量传递给与之关联的新分区；否则进入步骤6；若重新RCB分区，则步骤6-9所述分区即为新分区，否则即为老分区；
步骤6，根据步骤4更新的坐标，各进程对穿越分区的块体单元进行信息交换；
步骤7，各进程对分区内的块体单元进行接触检测，得到接触链表，与此同时计算当前时间步对应分区内且包含外重叠区每个块体单元所受的接触力和力矩；
步骤8，将各分区内重叠区中块体单元接触力和对应相邻分区的外重叠区中块体单元接触力进行交换；
步骤9，由刚体运动方程，根据接触力计算出每个块体单元当前时间步的加速度，再由verlet算法计算出每个块体单元下一时间步的速度；
步骤10，重复步骤4至步骤9计算下一时间步，直至计算完所有时间步，得到放矿后岩石散落的形态。
2.根据权利要求1所述基于并行RCB三维势函数离散元的放矿仿真方法，其特征在于，所述步骤2的具体过程如下：
2.1计算各块体单元的形心坐标，求块体单元各角点与形心之间的距离，取最大距离作为块体单元的半径，取半径的2.01倍作为块体单元的直径，将所有块体单元对应的直径中最大直径的1.05倍作为内、外重叠区的宽度；
2.2根据块体单元的编号，将块体单元平均分配到各进程中；
2.3创建保存初始分区信息的自定义类型cpuRCB_type，类型中包含四个变量，分别为：分区单元数量ne、单元索引列表liste、分区左下角坐标Domleftlow、分区尺寸Domsize；
2.3.1cpuRCB_type类型数组长度的确定
指定每个分区方向的分区层数，由RCB二分的性质得到：
(1)若第一个分区方向的分区层数是N1，则该方向最后一层分区条带数
总分区数
(2)若第二个分区方向的分区层数是N2，则该方向最后一层分区条带数
总分区数
(3)若第二个分区方向的分区层数是N3，则该方向最后一层分区条带数
总分区数
令N＝N1+N2+N3，则最后一层总分区数ncpu＝M1×M2×M3＝2N-3；
所有层的总分区数ncpuRCB＝L1×L2×L3，以ncpuRCB为大小创建cpuRCB_type类型数组；
2.3.2变量ne、liste、Domleftlow和Domsize的确定
(1)设定分区单元数量允许容差ratioRCB和二分法最大迭代次数niterRCB，指定二分方向为idimn，idimn＝1,2,3分别表示二分方向为X,Y,Z；
(2)以指定idimn方向的空间区域两侧坐标为边界，以中间坐标为界限将空间区域划分为两个区域，当二分方向idimn＝1时：初始空间区域左下角顶点坐标为(X0,Y0,Z0)，空间区域在三个方向上的长度分别为[Lx,Ly,Lz]，则该方向以X1＝X0+Lx/2作为分割线，分割为两个区域：
区域1：左下角坐标为(X0,Y0,Z0)，区域尺寸为[X1-X0,Ly,Lz]；
区域2：左下角坐标为(X1,Y0,Z0)，区域尺寸为[X0+Lx-X1,Ly,Lz]；
当块体单元的形心(Cx,Cy,Cz)都在区域1内时，即：
X0＜Cx＜X1-X0
Y0＜Cy＜Ly
Z0＜Cz＜Lz
若上述三个条件都满足，则判断该块体单位在区域1内，同理判断单元是否在区域2内，由此统计两个区域内的分区单元数量ne和单元索引列表liste；
(3)各进程归约统计得到区域1和2块体单元的实际总数量，若区域1块体单元总数量大于区域2块体单元总数量，则令新的分割线为X2＝(X0+X1)/2，反之，则令新的分割线为X2＝(X0+X1+Lx)/2；
(4)重复步骤(1)-(2)，直到区域1和2块体单元总数量的差值与总块体单元数量的比值小于允许容差ratioRCB时，或者循环次数大于最大迭代次数niterRCB时，退出循环，此时得到的分割线坐标Xn即为二分法要得到的值，且二分法后两个区域的空间位置：
区域1：左下角坐标为(X0,Y0,Z0)，区域尺寸为[Xn-X0,Ly,Lz]；
区域2：左下角坐标为(Xn,Y0,Z0)，区域尺寸为[X0+Lx-Xn,Ly,Lz]；
(5)将区域1的左下角坐标、尺寸分别存入cpuRCB_type类型数组的第1个元素中，区域2的左下角坐标、尺寸分别存入cpuRCB_type类型数组的第2个元素中；
(6)将上述(2)-(4)的操作称为idimn方向的第1层分区，则第2层分区为：将第1层分区得到的两个区域分别作为(2)中的初始空间区域，重复(2)-(3)，得到的4个子分区，存入cpuRCB_type类型数组的第3至第6个元素中；以此类推，直至idimn方向执行完第Nidimn层分区；
(7)按照上述(6)中的方法进行其他两个方向的分区，并填充到cpuRCB_type类型数组中；
2.4建立分区关联列表
2.4.1获取分区三维编号
分区三维编号由其所在空间分区拓扑中的位置确定，分区编号(I,J,K)代表该分区在所有分区中的排序为：第I行、第J列、第K层，其中，行代表X方向，列代表Y方向，层代表Z方向；
2.4.2确定进程编号和分区编号的关系
通过MPI_CART_CREATE函数创建笛卡尔拓扑通信域Cart_Comm，笛卡尔拓扑的各方向进程数量与分区在各方向的数量一致，由此，在通信域Cart_Comm中，由MPI_CART_RANK函数获取分区三维编号(I,J,K)所对应的进程编号，通过MPI_CART_COORDS函数由进程编号获取分区三维编号(I,J,K)；
2.4.3获取各分区的内外重叠区
将2.3中划分好的各分区的边界E向内收缩一个重叠区宽度，得到边界E内，将边界E向外扩展一个重叠区宽度得到边界E外，则边界E外到边界E之间的区域称为该分区的外重叠区，边界E到E内之间的区域称为该分区的内重叠区，边界E内以内的区域称为该分区的内区；
根据各分区中块体单元所在位置将块体单元标记为如下三种情况：
(1)若块体单元位于内区，标记为1；
(2)若块体单元位于内重叠区，标记为2；
(3)若块体单元位于外重叠区，标记为3；
2.4.4分区关联判断具体步骤：
(1)创建cpuDB_type类型的数组，数量为最后一层总分区数ncpu，cpuDB_type类型包含分区左下角坐标、分区尺寸、关联分区的个数ndlinkDB、关联分区列表dlinkDB；将2.3中cpuRCB_type类型数组第2N-1到2N-1个元素中分区的几何信息分别复制到cpuDB_type类型数组的第1到2N-1个元素中；
(2)确定分区三维编号到cpuDB_type类型数组的映射关系
分区的三维编号(I,J,K)与cpuDB_type类型数组中的编号icpu关系如下：
第I个面、第J个条带的第K个分区编号icpu是：
icpu＝(I-1)×M2×M3+(J-1)×M3+K
分区数组编号为icpu的分区，其三维编号(I,J,K)是：
I＝(icpu-1)/(M2×M3)+1
J＝(icpu-(I-1)×M2×M3-1)/M3+1
K＝icpu-(I-1)×M2×M3-(J-1)×M3
(3)筛选出有可能关联的分区
对三维编号为(I,J,K)的分区来说，有可能关联的分区为如下两类：
①一定会通讯的是：(I,J,K+1)
②可能会通讯的是：(I,J+1,1:M3)，(I+1,1:M2,1:M3)
其中，1:M3表示从1到M3的数组，1:M2表示从1到M2的数组；
(4)对可能关联的分区进行关联判断
对于cpuDB_type类型数组中编号为icpu的分区而言，要对所有可能关联的分区进行判别，具体判别方法如下：
①将可能通讯的分区三维编号(Ij,Jj,Kj)转化为cpuDB_type数组中编号jcpu；
②获取cpuDB_type数组中编号为icpu和jcpu的元素中的分区几何信息，即左下角坐标与尺寸，将左下角坐标的三个值分别加上三个方向的尺寸，求得分区三个方向的范围；对于分区序列(I,J+1,1:M3)中的分区，判断第三个方向是否在(I,J,K)分区第三个方向范围内，若在则分区关联，否则不关联；对于分区序列(I+1,1:M2,1:M3)中的分区，判断第二和第三个方向是否在(I,J,K)分区第二和第三个方向范围内，只要有一个方向不在范围内，则分区不关联；
2.5优化通讯批次
采用图论的方法优化通讯的顺序批次，保证每个批次中没有重复的通讯进程，具体步骤如下：
(1)标记cpuDB_type数组每个元素中所有关联分区为未使用状态，再将所有分区标记为未循环状态；
(2)从cpuDB_type数组中的第一个元素开始索引，如果该分区为未循环状态，则寻找该元素关联分区列表中第一个处于未使用且未循环状态的分区序号，如果有，则将其序号jcpu与该元素序号icpu记为一个通信对，最后将jcpu标记为循环且已使用状态；重复该操作，直到循环完所有元素，则完成了一批通讯批次；
(3)重复第(2)步，直到在某个批次中找不到通信对为止；
(4)按照求得的通讯批次的先后顺序，将cpuDB_type数组每个元素对应的关联分区编号重新填写；
(5)每一轮通信时，编号小的进程先接受后发送，编号大的进程先发送后接受。
3.根据权利要求1所述基于并行RCB三维势函数离散元的放矿仿真方法，其特征在于，步骤3所述时间步长小于等于允许的最大时间步长，且允许的最大时间步长为：
其中，Δt是允许的最大时间步长，ξ是阻尼比，
c是阻尼系数，m是块体单元的质量，k是刚度系数。
4.根据权利要求1所述基于并行RCB三维势函数离散元的放矿仿真方法，其特征在于，所述步骤4的具体过程如下：
当前时间步块体单元的速度，预估公式为：
更新坐标由下式进行：
其中，
为当前时间步预估速度，
为上一时间步速度，Δt为时间步长，
为上一时间步加速度，un+1为当前时间步位移，un为上一时间步位移。
5.根据权利要求2所述基于并行RCB三维势函数离散元的放矿仿真方法，其特征在于，所述步骤5的具体过程如下：
5.1判断是否需要进行重新RCB分区的标准如下：
(1)根据总的块体单元数/总分区数ncpu，得到平均单元数；
(2)对当前分区，若块体单元数小于平均单元数的某个比例，则认为该分区是坏的，赋值标记为1，否则为0；
(3)归约统计坏分区的数量；
(4)如果坏分区的数量大于总分区数ncpu的某个比例，则需要重新进行RCB分区；
5.2重新RCB方法：
5.2.1重复步骤2.3-2.4完成对空间区域的重新划分，得到新分区；
5.2.2新老分区关联
(1)判断是否关联
将每一个老分区与每一个新分区在三个方向比较其范围，若老分区在某个方向坐标的最小值大于新分区对应方向坐标的最大值，或在某个方向坐标的最大值小于新分区对应方向坐标的最小值，则说明该老分区与新分区无关联，否则有关联；
(2)将每个老分区关联的新分区编号记录到新老分区关联列表中；
5.2.3优化新老分区间的通讯批次
(1)分离新老分区关联列表：
对icpu号老分区的关联列表进行索引，若关联列表中有jcpu号新分区，再索引jcpu号新分区的关联列表，若存在icpu号老分区，则在jcpu号老分区的关联列表中剔除icpu号，并在新建的关联列表中jcpu号老分区关联的队列中添加-icpu号；
(2)将步骤2得到的两个关联列表调整为只包含比自身编号大的分区的状态：
icpu号老分区将其关联列表中的所有的分区标记为关联状态，再索引比icpu号大的老分区关联列表中是否有icpu号新分区，若有则将icpu号老分区标记关联状态，并将比icpu号大的老分区关联列表中的icpu号剔除，最后将所有老分区中处于关联状态的进程从小到大写入icpu号老分区关联列表中；由此方法按编号从小到大依次重填关联列表，即可得到关联列表只包含比自身编号大的分区的列表；
(3)重复步骤2.5，将上一步中所得的两个新关联列表分别优化，得到两个通讯批次；
(4)合并两个新关联列表，得到合并后的关联列表；
(5)按照5.2.2得到的新老分区关联列表标记合并后的关联列表中的发送和接受状态；
5.2.4新老分区单元信息互换
(1)创建dlinkDB_type类型数组dlinkDB，该数组长度为本老分区关联分区的数量，类型包含5个变量：关联的分区编号jcpu、需要传送给jcpu的单元索引列表listout、传送的单元个数neout、穿越到jcpu的单元列表listpass、穿越到jcpu的单元数量nepass；
(2)统计各分区内的块体单元数量及单元索引列表，分别存储到dlinkDB数组中的变量neout和listout中；
(3)按照5.2.3中合并后的通讯批次，按照dlinkDB中neout发送单元数量；
(4)按照上述(3)中的通讯批次，交换listout数组；
(5)按照索引列表listout将需要发送的单元的形心、速度、角速度打包到一个数组A中；
(6)将数组A发送到分区jcpu中，分区jcpu接受后，将接受到的数组解包，并填充到对应单元的形心、速度、角速度列表中。
6.根据权利要求1所述基于并行RCB三维势函数离散元的放矿仿真方法，其特征在于，所述步骤9的具体过程如下：
计算加速度的公式如下：
计算下一时间步速度的公式如下：
其中，
为当前时间步加速度，Fn+1为当前时间步接触力，m为块体单元质量，
为下一时间步速度，
为当前时间步预估速度，Δt为时间步长。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910397522.2A CN110211235B (zh) | 2019-05-14 | 2019-05-14 | 基于并行rcb三维势函数离散元的放矿仿真方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910397522.2A CN110211235B (zh) | 2019-05-14 | 2019-05-14 | 基于并行rcb三维势函数离散元的放矿仿真方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110211235A CN110211235A (zh) | 2019-09-06 |
| CN110211235B true CN110211235B (zh) | 2022-08-19 |
Family
ID=67787106
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910397522.2A Active CN110211235B (zh) | 2019-05-14 | 2019-05-14 | 基于并行rcb三维势函数离散元的放矿仿真方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110211235B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110780842A (zh) * | 2019-10-25 | 2020-02-11 | 无锡恒鼎超级计算中心有限公司 | 基于神威架构的船舶三维声弹性模拟计算的并行优化方法 |
| CN112052587B (zh) * | 2020-09-02 | 2023-12-01 | 中国人民解放军陆军工程大学 | 砂土垫层的三维细观离散体模型密实方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105912852A (zh) * | 2016-04-08 | 2016-08-31 | 河海大学 | 一种基于距离势函数任意凸多边形块体离散单元法 |
| CN106503365A (zh) * | 2016-11-03 | 2017-03-15 | 英特工程仿真技术（大连）有限公司 | 一种用于sph算法的分区搜索方法 |
| CN106528989A (zh) * | 2016-11-03 | 2017-03-22 | 英特工程仿真技术（大连）有限公司 | 一种分布式并行sph仿真方法 |
| CN106529146A (zh) * | 2016-11-03 | 2017-03-22 | 河海大学 | 一种基于距离势函数三维任意凸多边形块体离散单元法 |
-
2019
- 2019-05-14 CN CN201910397522.2A patent/CN110211235B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105912852A (zh) * | 2016-04-08 | 2016-08-31 | 河海大学 | 一种基于距离势函数任意凸多边形块体离散单元法 |
| CN106503365A (zh) * | 2016-11-03 | 2017-03-15 | 英特工程仿真技术（大连）有限公司 | 一种用于sph算法的分区搜索方法 |
| CN106528989A (zh) * | 2016-11-03 | 2017-03-22 | 英特工程仿真技术（大连）有限公司 | 一种分布式并行sph仿真方法 |
| CN106529146A (zh) * | 2016-11-03 | 2017-03-22 | 河海大学 | 一种基于距离势函数三维任意凸多边形块体离散单元法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110211235A (zh) | 2019-09-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10007742B2 (en) | Particle flow simulation system and method | |
| US5682322A (en) | Optimization processing for integrated circuit physical design automation system using chaotic fitness improvement method | |
| US8751556B2 (en) | Processor for large graph algorithm computations and matrix operations | |
| CN110211235B (zh) | 基于并行rcb三维势函数离散元的放矿仿真方法 | |
| CN113689556A (zh) | 一种块自适应型笛卡尔网格快速图映射方法及系统 | |
| CN115016951B (zh) | 流场数值模拟方法、装置、计算机设备和存储介质 | |
| US8554527B2 (en) | Particle simulator and method of simulating particles | |
| CN109522382A (zh) | 空间数据网格化统计方法及装置 | |
| CN113392587B (zh) | 一种大区域滑坡危险性评价的并行支持向量机分类方法 | |
| CN104459776A (zh) | 一种断裂分形特征优化计算方法 | |
| CN109917754B (zh) | 一种基于多种群分布估计算法的机器人装配单元多目标布局优化方法 | |
| US20230098447A1 (en) | Method for process allocation on multicore systems | |
| CN111414961A (zh) | 一种基于任务并行的细粒度分布式深度森林训练方法 | |
| KR101642823B1 (ko) | 이웃 탐색 연산 시스템 | |
| Cintra et al. | A parallel DEM approach with memory access optimization using HSFC | |
| JP5467262B2 (ja) | 粒子シミュレーション装置及び粒子シミュレーション方法 | |
| CN105787020A (zh) | 图数据划分方法及装置 | |
| CN115344383A (zh) | 一种基于进程并行的流线可视化并行加速方法 | |
| CN112989683A (zh) | 一种sph的向量化并行计算方法及装置 | |
| Nesmachnow et al. | Large-scale multithreading self-gravity simulations for astronomical agglomerates | |
| CN104407840B (zh) | 网格分解方法及系统 | |
| CN113760955B (zh) | 一种考虑盒子分形维数和产状的节理多因素分组方法 | |
| US20110119292A1 (en) | Nearest neighbor search method | |
| Tzovas et al. | Distributing sparse matrix/graph applications in heterogeneous clusters-an experimental study | |
| CN116088577B (zh) | 一种无人集群自主探索方法、系统、电子设备及介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |