CN106446433A - 一种求解材料大变形的双层混合并行最优输运无网格方法 - Google Patents

Abstract

本发明涉及一种求解材料大变形的双层混合并行最优输运无网格方法(parallel Optimal Transportation Meshfreep,pOTM)。pOTM是一种基于分布式多进程并行化与共享内存多线程并行化的双层混合并行架构的大规模并行方法，在MPI层发展了一套Shadow point/node计算数据划分方案，将OTM模拟进行分布式计算分解；在多线程层，发展了一套动态负载平衡的多线程并行化计算方案。pOTM方案有效的利用异构超级计算集群对OTM模拟进行线性甚至超线性加速，极大的提高了计算效率。

Description

一种求解材料大变形的双层混合并行最优输运无网格方法

技术领域

本发明涉及高效稳定的求解包含极大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型、多相变及与加载历史相关的材料动态响应的连续介质力学问题的数值方法，特别是涉及一般连续介质力学问题的无网格求解方法，它是一种结合了分布式多进程并行化与共享内存多线程并行化的双层混合并行最优输运无网格计算方法，属于计算力学领域。

背景技术

基于拉格朗日网格和欧拉网格的传统数值方法，如有限单元法、有限体积法及有限差分法等，因其本身固有的网格特性限制了它们在极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等愈加复杂的实际工程问题方面的应用。对于基于拉格朗日网格的求解方法，由于网格的存在所引起的困难主要体现在：①网格大变形引起收敛缓慢、结果不准确，甚至过早的终止计算；②网格尺寸过小时，时间步趋于零，将导致计算过程崩溃；③网格重构的困难及进行新旧网格之间数据插值计算而引入的扩散误差和错误。对于基于欧拉网格的求解方法，由于网格的存在所引起的困难主要体现在：①网格固定材料在运动，分析中难以处理与加载历史相关的材料动态响应以及裂纹扩展问题；②欧拉方法仅对经过网格单元边界的质量、动量及能通量进行跟踪，使得自由表面的位置，可变形边界和移动材料界面难以确定；③复杂的网格生成程序，将问题域不规则的几何形状转换到规则的计算域中，有时会带来昂贵的数值映射开销。

基于网格方法的所有困难和局限性皆源自于网格的存在，为了解决这个困难，下一代数值计算方法—无网格方法得到了广泛快速的发展，无网格方法的目标是通过以节点的方式对整个问题域进行离散，从而完全(至少从部分上)消除对网格的依赖性。其主要的优越性体现在：①计算过程中不涉及网格，不存在网格畸变，网格重构等困难，非常适合处理大变形问题；②H型和P型自适应实现变得相对简单；③可以很容易的构造任意阶次连续性的形函数；④收敛结果往往大大优于基于网格的方法；⑤体积自锁问题可以通过调整形函数扩张参数得到缓解。然而，无网格方法的发展也正面临着亟待解决的困难：①绝大部分无网格插值函数都不满足克罗内克属性，使得施加位移边界条件变得非常困难；②在缺少网格的条件下，在对无网格伽辽金法进行“弱”形式下的等效积分时会有较大困难，大量研究文献表明，几种可能的积分方案(节点积分方案、应力点积分方案、背景网格积分方案)都不同程度的面临着处理拉应力不稳定性、数值积分误差等困难；③缺乏严密的收敛性与误差理论的数学分析；④由于无网格方法中使用高阶插值函数，对高性能计算提出了很大的挑战。

新兴的最优输运无网格方法(Optimal Transportation Meshfree method，OTM方法)通过有机的结合物质点空间离散、局部最大熵插值函数和最优输运理论时间离散，有效克服了基于网格和无网格方法遇到的诸多困难，其优越性主要体现在：①采用物质点和节点对问题域进行离散，物质点充当数值积分点，避免了计算结果在拉应力不稳定性；②采用局部最大熵插值函数满足“弱”克罗内克属性，使得位移边界条件难施加的困难得以解决；③最优运输理论保证了其时间离散形式在理论上严格满足质量与动量守恒定律，具有严格的收敛性数学证明；④动态调整物质点与节点的连接关系自动进行动态接触面的检测以及接触力的计算，没有额外的计算成本。OTM方法的优异特性使其在求解包含超大变形、任意边界条件、复杂几何结构，与加载历史相关的材料动态响应的连续介质力学问题方面有巨大的潜力。

OTM方法是一套显式增量更新的拉格朗日无网格方法，与显式有限单元求解法一样，在求解高非线性动力学问题时，OTM方法同样面临着亟需解决的大计算量问题。在OTM方法中采用物质点与节点相结合的方式对问题域进行离散，其中计算邻域的运动学信息，包括位移、速度、加速度与温度等被存储在节点上，而材料的物理信息，比如变形、应力、材料内部参数等存储于物质点上。在求解过程中主要的计算负荷集中于预估物质点上的材料响应，即节点力的求解；同时在求解时间历程中随着物质点邻域的动态更新，物质点形函数及其导数重新计算是另一耗时的过程。因此如何提高OTM方法在解决工程系统中极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题时的计算效率是OTM方法需要重点解决的问题之一，由于在OTM方法中主要的计算都是在物质点上进行的，各个物质点之间没有直接的数据交换，非常适宜对其计算过程进行并行化处理，通过成功运用并行计算机及并行策略，实现节点力、形函数及导数计算过程的并行化，有望为OTM方法的求解过程节省可观的计算开销，大幅提升其计算效率。

发明内容

本发明技术解决方案：克服现用技术的不足，提出了一种求解材料大变形的双层混合并行最优输运无网格方法(parallel Optimal Transportation Meshfreep,pOTM)，pOTM方法结合了分布式多进程并行化与共享内存多线程并行化，有效的利用异构超级计算集群对OTM模拟进行线性甚至超线性加速。分为分布式多进程并行化和多线程动态负载平衡，其中：

分布式多进程并行化步骤如下：

设Ω表示d维的连续介质问题域，被离散为一组物质点集{x_p,k,p＝1,2,…,M；k＝0,1,…,n}和一组节点集{x_a,k,a＝1,2,…,N；k＝0,1,…,n}，采用pOTM方法(parallelOptimal Transportation Meshfreep,pOTM,双层混合并行最优输运无网格方法)对连续介质域内的材料动态响应进行求解时，

第一步，定义总计算步数n，定义分布式处理器数量在t_k＝0时刻根据处理器数量P，对连续介质问题域中的物质点进行划分并将划分后的物质点集分别发送到各个处理器中；

第二步，计算t_k时刻处理器上各个物质点的邻域，得到各个物质点的邻域范围

第三步，根据第二步各物质点邻域的边界范围，确定处理器的节点RangeBox，由Range Box得到该处理器物质点邻域内的所有节点，

第四步，根据第三步确定的Range Box，将各个Range Box之间相互重叠的部分定义为Shadow Box，处于Shadow Box的节点被定义为Shadow node，得到的Shadow node，将它们存储至共享数据交换表C_k；

第五步，初始化t_k时刻物质点形函数N_a,k(x_p,k)、物质点形函数导数物质点变形梯度F_p,k；

第六步，根据第一步至第五步，完成了进程层各个处理器的数据分配工作，接下来将转入每个单独的处理器在每个处理器中采用多线程动态负载平衡计算方案，得到材料的动态响应；

多线程动态负载平衡步骤如下：

第七步，根据第六步，计算过程从分布式多进程层转入每个单独的处理器进行多线程并行计算，在此步骤中用户首先定义线程数量，然后各处理器需要检查实际可用的线程数是否满足所需求的线程数，不满足则按实际可用的最大线程数为最终调用的线程数；

第八步，根据第七步确定的线程数，将处理器上所包含的物质点集M^I动态负载平衡分配到各个线程

第九步，各处理器上各线程根据第八步分配得到的物质点子集进行求解计算，得到各处理器的局部质量矩阵局部节点力及局部加速度

第十步，各个处理器将所包含的共享节点的局部节点力局部质量矩阵局部加速度发送至数据交换表C_k中，并接收该共享节点在其他各个处理器计算结果，完成数据同步，进行数据组装得到全局质量矩阵全局节点力和全局加速度

第十一步，更新节点数据：利用第十步的结果进行数据更新，得到t_k+1时刻的节点坐标节点速度

第十二步，t_k+1时刻，更新物质点坐标更新物质点变形梯度F_p,k+1、更新物质点本构、更新物质点邻域更新物质点形函数N_a,k+1(x_p,k+1)、更新物质点形函数导数

第十三步，根据第十二步更新的物质点邻域更新各Range Box大小，搜索并更新Range Box内包含的节点集，

第十四步，根据第十三步确定t_k+1时刻各个处理器的Range Box范围，重新搜索Shadow node并更新共享数据交换表C_k；

第十五步，判断时间步t_k+1，若t_k+1＝t_n，代表已计算至最后一步，此时退出计算，最终得到物质点和节点从t₀→t_n时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析；若t_k+1≠t_n，则转入第七步进行下一次迭代计算，直至时间步t_k+1＝t_n为止，最终得到物质点和节点从t₀→t_n时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析。

本发明与现有技术相比的优点在于：

1)本发明所提的pOTM方法是OTM方法的并行化方法，在pOTM中将物质点和节点划分成子集分配到不同的计算服务器中，通过使用信息传递接口(Message PassingInterface,MPI)进行不同服务器中的数据交流，实现了对OTM计算方法的并行化，有效的利用异构超级计算集群对OTM模拟进行线性甚至超线性加速。

2)在pOTM方法中，提出了一种Shadow point/node计算数据划分方案：根据处理器分配得到的物质点进行物质点邻域计算，由上所有物质点邻域的边界组成一个邻域范围，这个邻域范围在本发明中被定义为Range Box。Range Box里包含了上所有物质点邻域范围内的节点，这个Range Box可以理解为原问题域的一个子域，基于这个Range Box每个处理器就可以单独完成子域的计算。而不同Range Box之间的会有相互重叠的部分，将各个Range Box之间相互重叠的部分定义为Shadow Box，处于Shadow Box的物质点/节点被定义为Shadow point/node，其中Shadow nodes代表各个Range Box之间的共享节点，将它们存储至共享数据交换表C_k。每个处理器完成单独的计算之后，通过MPI对C_k中共享节点数据进行数据交互，最后得到原问题域的求解结果。

3)在pOTM方法中，对于每个子域在处理器上的计算，提出了一种动态负载平衡的多线程并行化方案：在t_k时刻将处理器上的物质点均分给各线程，完成一次计算后得到各个物质点的计算量大小，而在下一个时间步t_k+1，则根据t_k时刻预估的计算量对各个线程上物质点按照计算量的大小重新进行分配，以此来保证每个线程获得的计算量所需要的计算时间相当，实现了每步迭代计算中子域的动态分解，避免一般随机分配计算量的做法所引起的负载不均问题。同时，多线程动态负载平衡算法不涉及任何数据交换与在内存空间中的传递，不会引发任何附加计算量。

附图说明

图1为本发明一种连续介质力学问题中求解材料动态响应的双层混合并行最优输运无网格方法流程图；

图2为本发明通过物质点与节点相结合的方式对几何模型进行空间离散的示意图；

图3为本发明Range Box的定义方案，图示将问题域中的物质点划分成3份分配给3个处理器，各处理器确定各自的Range Box，各Range Box重叠部分为Shadow Box；

图4为本发明多线程动态负载平衡方案，图示中在MPI层将物质点划分成两份，分别分配到两个处理器P1和P2中进行处理，在线程层(PThreads)中，又将物质点划分成4份，分别分配给4个线程(T1、T2、T3、T4)进行计算；

图5为本发明在t_k→t_k+1，计算更新节点运动学信息之后MPI层数据交互方案；

图6为本发明的仿真算例，右图是采用pOTM方法模拟Al6061-T6金属板被440c钢球以2.7km/s的速度撞击效果图，该仿真模拟包含1,116,024个物质点，683,772个自由度，在1024个计算核心的LLNL’s Hera服务器上进行计算；左图表示每个处理器上分别采用单线程、2线程、4线程、8线程的加速效果。

具体实施方式

下面结合附图，对本发明一种连续介质力学问题中求解材料动态响应的双层混合并行最优输运无网格方法的技术方案做进一步说明。

如图1所示，本发明实现如下：

分布式多进程并行化步骤如下：

设Ω表示d维的连续介质问题域，被离散为一组物质点集{x_p,k,p＝1,2,…,M；k＝0,1,…,n}和一组节点集{x_a,k,a＝1,2,…,N；k＝0,1,…,n}，采用pOTM方法(parallelOptimal Transportation Meshfreep,pOTM,双层混合并行最优输运无网格方法)对连续介质域内的材料动态响应进行求解时，

第一步，定义总计算步数n，定义分布式处理器数量在t_k＝0时刻根据处理器数量P，对连续介质问题域中的物质点进行划分并将划分后的物质点集分别发送到各个处理器中；

在OTM方法中采用物质点与节点相结合的方式对所求问题域进行离散，物质点和节点的生成方式及初始位置可由用户根据不同的算法来确定，比如随机插入物质点与节点，以距离物质点最近的n+1个节点作为其计算邻域(n为问题域的维度)，计算邻域的运动学信息，其中包括位移、速度、加速度与温度等，将被存储在节点上。材料的物理信息，比如变形、应力、材料内部参数等，将存储在物质点上。在本发明的具体实现中，为了保证计算的精度，采用三角形单元(二维情况)或者四边形单元(三维情况)对节点进行网格划分，每个单元的形心取为物质点。如图2所示，空心点代表节点：x_a,k(a为离散域中节点索引号，代表第几个节点，k代表第几个时间步，时间步长及步数由用户控制)，实心点代表物质点：x_p,k(p为离散域中物质点索引号，代表第几个物质点，k代表第几个时间步)。

通过上述方式获得了问题域离散化的节点集与物质点集，在本步骤中即是要对物质点集进行划分，划分成不同的子集，并将不同的子集分配到不同的处理器进行分布式计算，本发明采用Metis算法对物质点集划分，后续计算在物质点上进行不涉及网格。

第二步，计算t_k时刻处理器上各个物质点的邻域，得到各个物质点的邻域范围其中d_p,k为物质点的动态邻域大小，在本发明中由公式(1)确定：

d_p,k＝δx×h (1)

其中δx为人工系数，取值范围为2.0～5.0，h为物质点的大小，物质点的大小通过公式(2)确定：

其中d为所求问题域的维度，V_t为物质点的体积，在步骤2)通过初始网格可以得到物质点的初始体积，随后网格被抛弃，物质点的体积通过在每个时间步求解质量守恒方程动态计算而得，由此计算出物质点的大小h。

第三步，根据第二步各物质点邻域的边界范围，确定处理器的节点RangeBox，由Range Box得到该处理器物质点邻域内的所有节点，

具体为：对处理器上的各物质点，计算和 其中L_k代表处理器上物质点邻域的下边界，U_k代表处理器上物质点邻域的上边界，定义以L_k和U_k为对角顶点的长方形(二维情况)或者长方体(三维情况)为Range Box，在Range Box中包含了该处理器物质点邻域内的所有节点， 如图3所示有3个MPI处理器，将问题域中的物质点划分成3份分配给每个处理器之后，由其上的物质点邻域分布确定了各自的Range Box，位于Range Box范围内的节点即为处理器t_k时刻的节点集，

第四步，根据第三步确定的Range Box，将各个Range Box之间相互重叠的部分定义为Shadow Box，处于Shadow Box的节点被定义为Shadow node，得到的Shadow node，将它们存储至共享数据交换表C_k；

第五步，初始化t_k时刻物质点形函数N_a,k(x_p,k)、物质点形函数导数物质点变形梯度F_p,k；

第六步，根据第一步至第五步，完成了进程层各个处理器的数据分配工作，接下来将转入每个单独的处理器在每个处理器中采用多线程动态负载平衡计算方案，得到材料的动态响应；

多线程动态负载平衡步骤如下：

第七步，根据第六步，计算过程从分布式多进程层转入每个单独的处理器进行多线程并行计算，在此步骤中用户首先定义线程数量，然后各处理器需要检查实际可用的线程数是否满足所需求的线程数，不满足则按实际可用的最大线程数为最终调用的线程数；

第八步，根据第七步确定的线程数，将处理器上所包含的物质点集M^I动态负载平衡分配到各个线程

在本发明中，为线程层的物质点划分发展了一套动态负载平衡方案，具体为：假设可用的线程数为t_k时刻每个线程分配到的物质点为其中代表处理器上的总物质点，令C(x_p,k)代表物质点x_p,k在t_k时刻的计算量，于是分配给线程的计算量可以表达为公式(3)所示：

于是处理器在t_k时刻的总计算量为：

由于从t_k→t_k+1，物质点只有非常微小的变化，因此在t_k+1时刻，各线程的物质点可以通过处理器在t_k时刻的计算量进行重构，即：

如图4所示，首先是在MPI层将物质点划分成两份，分别分配到两个处理器P1和P2中进行处理。然后在线程层(PThreads)中，又将物质点划分成4份，分别分配给4个线程(T1、T2、T3、T4)进行计算。通过动态负载平衡多线程并行化方案，实现了每步迭代计算中子域的动态分解，避免一般随机分配计算量的做法所引起的负载不均问题，如图6所示为本发明的仿真算例，采用pOTM方法对440c球形弹丸以2.7km/s的速度撞击Al6061-T6金属板的进行仿真模拟，该问题由1,116,024物质点及683,772个自由度，通过1024核心服务器进行计算，对比了在MPI层上，不同线程数(单线程、2线程、4线程及8线程)对计算加速因子的影响；

第九步，各处理器上各线程根据第八步分配得到的物质点子集进行求解计算，得到各处理器的局部质量矩阵局部节点力及局部加速度

第十步，各个处理器将所包含的共享节点的局部节点力局部质量矩阵局部加速度发送至数据交换表C_k中，并接收该共享节点在其他各个处理器计算结果，完成数据同步，进行数据组装得到全局质量矩阵全局节点力和全局加速度其中Q∈C_k即为ShadowBoxes中的共享节点。如图5所示，在t_k→t_k+1，计算更新节点运动学信息之后，对RangeBox中的节点进行判断，如果该节点属于Shadow node则将数据信息发送到数据交换表C_k，同时从C_k接收来自其他处理器的计算结果，进行节点运动学数据的组装，C_k也在每个时间步得到更新；

第十一步，更新节点数据：利用第十步的结果进行数据更新，得到t_k+1时刻的节点坐标节点速度

第十二步，t_k+1时刻，更新物质点坐标更新物质点变形梯度更新物质点本构关系、更新物质点邻域更新物质点形函数更新物质点形函数导数

第十三步，根据第十二步更新的每个物质点的邻域搜集处理器上各个物质点邻域内的节点，将它们更新为处理器t_k+1时刻的节点集，

第十四步，根据第十三步确定了t_k+1时刻物质点的邻域，定义Range Box的边界坐标值L_k+1和U_k+1得到更新，因此Range Box在每个时间步的迭代计算中将会动态更新，使得共享节点集C_k中的共享节点随着计算的进行将会被动态更新；

第十五步，判断时间步t_k+1，若t_k+1＝t_n，代表已计算至最后一步，此时退出计算，最终得到物质点和节点从t₀→t_n时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析；若t_k+1≠t_n，则转入第七步进行下一次迭代计算，直至时间步t_k+1＝t_n为止，最终得到物质点和节点从t₀→t_n时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (3)

1.一种求解材料大变形的双层混合并行最优输运无网格方法，其特征在于：分为分布式多进程并行化和多线程动态负载平衡，其中：

分布式多进程并行化步骤如下：

设Ω表示d维的连续介质问题域，被离散为一组物质点集{x_p,k,p＝1,2,…,M；k＝0,1,…,n}和一组节点集{x_a,k,a＝1,2,…,N；k＝0,1,…,n}，采用pOTM方法(parallelOptimal Transportation Meshfreep,pOTM,双层混合并行最优输运无网格方法)对连续介质域内的材料大变形问题进行求解时，

第一步，定义总计算步数n，定义分布式处理器数量I＝1,…,P，在t_k＝0时刻根据处理器数量P，对连续介质问题域中的物质点进行划分并将划分后的物质点集分别发送到各个处理器中；

第二步，计算t_k时刻处理器上各个物质点的邻域，得到各个物质点的邻域范围

第三步，根据第二步各物质点邻域的边界范围，确定处理器的节点Range Box，由Range Box得到该处理器物质点邻域内的所有节点，

第四步，根据第三步确定的Range Box，将各个Range Box之间相互重叠的部分定义为Shadow Box，处于Shadow Box的节点被定义为Shadow node，得到的Shadow node，将它们存储至共享数据交换表C_k；

第五步，初始化t_k时刻物质点形函数N_a,k(x_p,k)、物质点形函数导数物质点变形梯度F_p,k；

第六步，根据第一步至第五步，完成了进程层各个处理器的数据分配工作，接下来将转入每个单独的处理器在每个处理器中采用多线程动态负载平衡计算方案，得到材料的动态响应；

多线程动态负载平衡步骤如下：

第七步，根据第六步，计算过程从分布式多进程层转入每个单独的处理器进行多线程并行计算，在此步骤中用户首先定义线程数量，然后各处理器需要检查实际可用的线程数是否满足所需求的线程数，不满足则按实际可用的最大线程数为最终调用的线程数；

第八步，根据第七步确定的线程数，将处理器上所包含的物质点集M^I动态负载平衡分配到各个线程

第九步，各处理器上各线程根据第八步分配得到的物质点子集进行求解计算，得到各处理器的局部质量矩阵局部节点力及局部加速度

第十步，各个处理器将所包含的共享节点的局部节点力局部质量矩阵局部加速度发送至数据交换表C_k中，并接收该共享节点在其他各个处理器计算结果，完成数据同步，进行数据组装得到全局质量矩阵全局节点力和全局加速度

第十一步，更新节点数据：利用第十步的结果进行数据更新，得到t_k+1时刻的节点坐标节点速度

第十二步，t_k+1时刻，更新物质点坐标更新物质点变形梯度F_p,k+1、更新物质点本构、更新物质点邻域更新物质点形函数N_a,k+1(x_p,k+1)、更新物质点形函数导数

第十三步，根据第十二步更新的物质点邻域更新各Range Box大小，搜索并更新Range Box内包含的节点集，

第十四步，根据第十三步确定t_k+1时刻各个处理器的Range Box范围，重新搜索Shadow node并更新共享数据交换表C_k；

第十五步，判断时间步t_k+1，若t_k+1＝t_n，代表已计算至最后一步，此时退出计算，最终得到物质点和节点从t₀→t_n时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析；若t_k+1≠t_n，则转入第七步进行下一次迭代计算，直至时间步t_k+1＝t_n为止，最终得到物质点和节点从t₀→t_n时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析。

2.根据权利要求1所述的一种求解材料大变形的双层混合并行最优输运无网格方法，其特征在于：所述第四步中共享数据交换表C_k的定义方法为：

(1)计算处理器上各个物质点的邻域得到各个物质点的邻域

(2)根据(1)计算所得的邻域确定处理器的节点Range Box，得到各个处理器的Range Box；

(3)根据(2)计算得到的各个处理器I＝1,…,P，的Range Box，将各个Range Box之间的重叠部分定义为共享节点Shadow node，得到共享节点Shadow node；

(4)根据(3)中所得的共享节点Shadow node，将共享节点Shadow node保存至共享数据交换表C_k，得到共享数据交换表C_k；基于此共享数据交换表C_k，在每个处理器完成单独的计算之后，通过MPI对C_k中共享节点数据进行分布式进程间的数据交互，最后得到原问题域的求解结果。

3.根据权利要求1所述的一种求解材料大变形的双层混合并行最优输运无网格方法，其特征在于：所述第八步将处理器上所包含的物质点集M^I动态负载平衡分配到各个线程实现步骤如下：

设可用的线程数为i＝1,2,…,N，t_k时刻每个线程分配到的物质点为其中代表处理器上的总物质点，令C(x_p,k)代表物质点x_p,k在t_k时刻的计算量，

(1)根据每个物质点x_p,k在t_k时刻的计算量C(x_p,k)，得到t_k时刻分配给线程的计算量

(2)根据(2)所得每个线程的计算量得到处理器在t_k时刻的总计算量

(3)根据(2)所得处理器在t_k时刻的总计算量预估处理器上每个线程在t_k+1时刻的计算量

(4)根据(2)所得对物质点进行重新分组得到通过满足来获得动态负载平衡多线程并行化方案。