基于OpenMP的近场动力学模拟并行实现

徐凤洲，张健飞

(河海大学 力学与材料学院，江苏 南京 211100)

0 引 言

材料和结构破坏的计算机数值模拟一直是传统数值方法面临的一个难题，其局限性源自连续介质力学理论中的连续性假设和实际问题不连续之间的矛盾，在分析裂纹问题时，裂纹的起裂和扩展都必须引入外部的准则进行判断，无疑增加了分析问题的复杂性。近场动力学(peridynamics，PD)作为一种新兴的基于非局部作用的理论[1]和方法，在分析不连续问题时，有其独特的优势。近场动力学的运动(平衡)方程是采用空间积分形式，并基于牛顿第二定律，因此，对连续体和非连续体均适用。因其本构模型中已经包含了对损伤和断裂的描述，所以不需要附加开裂判据和裂纹路径分析，裂纹自然而然地萌生并扩展[2]。然而，近场动力学在分析裂纹破坏等不连续问题的过程中，需要将模型离散为一系列包含物质信息的物质点，进一步需要搜索每个物质点邻域范围内的，与其相互作用的其他物质点，同时，在求解计算过程中，为了保证数值准确性和稳定性还需要进行大量的迭代求解[3]。因此，近场动力学在实际数值模拟中面临着计算复杂度高、计算量大等问题[3-4]，这无疑会有碍于近场动力学的进一步发展和实际应用。

并行计算是目前扩大计算规模和加快计算速度的重要途经，开展PD数值模拟的并行计算研究可以有效扩大计算规模和提高计算速度。目前，PD数值模拟的并行计算研究采用的方法主要有MPI并行[5]、基于GPU和 CUDA的多线程并行[3]和基于GPU和OpenCL的多线程并行[6]等，但是这些方法的实现过程较为复杂，且对计算设备要求较高。OpenMP(open multi-processing)是一套用于共享内存并行系统多线程程序设计的指导性注释(compiler directive)，是为共享内存的多处理器系统设计的并行编程方法。OpenMP特别适合用在多核处理器计算机上运行的并行程序设计，可以快速实现PD模拟的并行化。目前，基于OpenMP的PD并行计算研究较少，缺乏具体实现过程描述。本文针对PD计算特点，对PD程序结构做归纳和总结，并基于OpenMP对PD数值模拟的并行实现方法进行研究，对程序的主要模块做并行化，编制并行程序，通过经典算例测试程序的性能，以期能提高计算效率。

1 基于键型的PD串行算法

近场动力学理论目前主要分为键型PD理论[7]和态型PD理论[8]。考虑到键型理论较早提出[7]，且研究比较深入和成熟[2]，本文选择键型PD理论，并采用键型PMB模型。PD理论的计算框架[9]主要包括模型离散、数据输入和存储、初始和边界条件施加、时间积分、物质点内力求解、破坏模型、数据输出等。另外，模型离散后需要根据近场尺寸范围搜寻并构建邻接节点域[9]。根据已有的PD计算框架，本文首先采用C++程序语言编制键型PD串行程序，并采用准静态模拟和动态裂纹扩展模拟2种模型。前者适用于不考虑损伤情况的一般弹性变形问题，后者则适用于快速加载条件下的裂纹扩展问题。以动态裂纹扩展模拟为例，PD程序的大体结构如图1所示。

图1 PD程序结构图

1.1 准静态模拟

对于准静态模拟，本文采用自适应动态松弛法求解问题，不考虑物质点的损伤情况，适用于模拟一般弹性变形过程，例如板的轴向受拉问题等。对于准静态问题，PD的运动(平衡)方程为

(1)

在准静态模拟中，程序按执行顺序大体可以分为以下几个模块：

(1)读取模型信息，将模型离散化——Data input模块。

(2)构建每个物质点的邻接节点域——Neighbor模块。

(3)计算每个物质点的表面修正系数，改善PD模型自身的边界效应问题——Compute kave模块。

(4)参数初始化——Initial模块。

(5)时间积分，更新物质点位置。其中又包含边界条件施加、计算物质点内力密度和自适应阻尼项、检查收敛性、更新物质点位移值，结果输出等子模块——Time integration模块。

1.2 动态裂纹扩展模拟

在模拟裂纹扩展问题中，需要考虑物质点的损伤情况，主要有2种类型：一种是准静态裂纹扩展；一种是动态裂纹扩展。本文采用动态裂纹扩展模拟，这是因为动态裂纹扩展应用较为广泛，与准静态模拟不同，它是采用显式求解方法，并且裂纹分叉结果更加明显、更加直观、可视化更好。对于动态问题，PD的运动(平衡)方程为

(2)

式中各参数的意义与准静态问题相同。在动态裂纹扩展模拟中，程序的大体结构与上述准静态模拟相似，但由于求解方法不同，在时间积分模块(Time integration)中，去掉了计算自适应阻尼项、检查收敛性等子模块，并且在计算物质点内力密度子模块中，增加了表征物质点损伤的量，需要分别对键破坏以及物质点损伤作出判断和计算。

2 OpenMP并行编程模型

OpenMP是一种面向共享内存以及分布式共享内存的、多处理器多线程并行语言。其编程模型是以线程为基础，通过编译制导语句在原来串行程序的基础上显示进行并行化处理[10]。对原有串行程序的改动较小，具有操作简便、开发效率高等优点。OpenMP采用的执行模式是串行→并行→串行→并行→串行……，如图2所示。这种执行模式的核心在于并行区域中线程的派生/缩并(Fork/Join)[11-12]。

图2 OpenMP并行程序的运行过程

3 基于OpenMP的PD并行算法

为了对已编制的PD串行程序进行更好的并行化处理，本文首先测试了PD串行程序各模块的计算时间，再对计算耗时长的模块进行并行化。测试结果显示,计算耗时长的模块主要集中在搜索每个物质点近场范围内的邻域节点(Neighbor模块)、计算每个物质点的表面修正系数(Compute kave模块)，以及时间积分(Time integration模块)这3个模块，并进一步考查了单个时间步内，时间积分模块(Time integration)中各子模块的计算时间，结果发现，计算物质点内力密度子模块的计算耗时占了这一模块(Time integration)计算耗时的大部分(80%左右，计算设备的性能不同，结果会存在差异)，所以在对程序进行并行处理时，需要着重对这一部分进行并行优化。分析了程序的热点循环和热点函数后，在对应的模块内，主要针对无循环依赖的循环体，增添OpenMP编译制导语句，对程序进行并行化。

3.1 对Neighbor模块的并行化

程序中Neighbor模块的主要作用是搜索每个物质点近场范围内的，与其相互作用的其他物质点，构建邻接节点域。在程序的这一模块中，将构建3个数组分别存放所有物质点的邻域节点号(按物质点排列顺序存放)(数组1)(Neighbor)、所有物质点的邻域节点个数(数组2)(numNeighbor)，以及所有物质点的起始邻域节点在数组1中的下标位置(数组3)(Nodefam)。由于数组1是存放所有物质点的邻域节点号，因此,必须按照物质点的排列顺序，依次将每个物质点的邻域节点号顺序存放到数组1中。数组3中存放的是所有物质点的起始邻域节点在数组1中的下标位置，由于数组1是按照物质点排列顺序存放的，因而数组3也需要按照顺序计算和存放数据。

整个Neighbor模块将由两个大的循环体构成。第一个循环体是二重嵌套循环，主要完成以下过程：搜索每个物质点的邻域内节点，构建邻接节点域；确定每个物质点的邻域内节点个数并存放到数组2中。由于第一个循环体没有循环依赖，每个物质点可以独立运行，因而可以进行并行化，并行化主要过程如图3所示。第二个循环体构建了上述的数组1和数组3，由于数组1和数组3中数据的存放需要按照物质点排列顺序，即循环指标的先后次序，排队执行，因而无法进行并行化，只能依照原有的串行方式执行。好在第二个循环体的运行时间相比第一个循环体要小得多，一般只有零点几秒，所以并没有影响到整个Neighbor模块的并行优化效果，这一模块的部分并行代码段如图4所示。private()括号内为声明的私有变量列表，需填上循环体内的私有变量(特别是进行写操作的变量)，确保每个物质点在程序执行过程中互不干扰，避免出现数据竞争现象。并行语句for指令后不带schedule子句，系统会默认采用静态调度方式(调度过程见图5)，即如果将循环迭代的总次数设为n，并行区域内线程总数为p，那么分配给每个线程约n/p次连续的迭代[13]。因此不同处理器在对共享数组(例如数组2)进行写操作时，操作的是数组中相隔比较远的元素，从而避免了伪共享[13]的问题。

图3 Neighbor模块并行化示意图(以3个线程为例)

图4 并行代码段1(构建物质点邻接节点域)

3.2 对Compute kave 模块的并行化

Compute kave模块的主要作用是改善PD模型自身的边界效应，对物质点的微模量系数值作出相应的修正。为了保证修正系数的准确性，将这一部分迭代了20次[14]。每一次迭代，每个物质点需要对邻域内节点做一次循环，所以是一个二重嵌套循环。OpenMP可以对嵌套循环体内的任意一个循环进行并行化，为了避免过多的并行开销导致并行效率降低，本文选择对外层循环做并行优化，而内层循环依然保持串行执行方式。在并行语句for指令后增添schedule(guided)子句，系统将会采用指导性调度方式进行调度(调度过程见图6)。指导性调度方式是一种动态方式，与静态调度方式(schedule(static))不同，分配给每个线程的迭代块不是均等的，而是在开始时迭代块会比较大，后来逐渐减小，从而有效地减少了调度开销[13]。这样处理是考虑到每个物质点的邻域内节点个数存在差异，特别是边界处的物质点，其邻域内节点个数会大大减少，如果仍然采用静态调度方式，则难以保证线程在程序执行过程中的负载平衡，而动态调度方式(schedule(dynamic))又存在额外的开销，不利于提高并行效率，所以最终考虑采用指导性调度方式(schedule(guided))。如果需要连续对多个for循环体进行并行化，那么在并行语句#pragma omp parallel下方需要增添一对大括号将整个并行区域括起来。

3.3 对Time integration模块的并行化

Time integration模块一般是整个程序中计算时间最长的部分(具体运行时间与设置的迭代步数有关)。这一模块的主体是一个大的迭代循环，可能包含成千上万次的迭代。迭代过程中，前后两次迭代是存在相互依赖关系的，并且迭代循环的次数在某些情况下也是无法提前确定的，所以对迭代循环整体无法运用OpenMP进行并行优化，但是在每一次迭代时间步下，程序的执行过程是可以并行化的。本文主要对Time integration模块中的计算物质点内力密度，更新物质点位移和速度等子模块进行并行优化，在程序的每一次迭代中这些子模块都会相应地被执行一次。计算物质点内力密度子模块的部分并行代码段如图7所示。对涉及物质点邻域节点循环的for循环体，本文都是采用指导性调度方式，而对于一般的物质点循环，例如更新物质点位移、速度等都是采用静态调度方式。

图7 并行代段码2(计算物质点内力密度)

在整个程序的并行化和调试过程中会用到一些OpenMP运行环境操作函数，例如omp_set_num_threads(n)，此函数用来确定并行区域内运行的子线程数量，只能在并行区域之外调用[13]。实际操作中，一般将此函数放在编译制导语句#pragma omp parallel的前面，通过设定n值确定并行区域运行时子线程的数量。

4 计算结果与分析

为了更好地考查程序的并行性能，以及计算规模与并行加速比之间的关系，选择以下3个算例，物质点总数分别为10 600，103 000，253 000。通过对比串行程序与并行程序的运行计算结果，验证并行程序的正确性，计算并分析3个算例对应的并行加速比情况。下面将详细阐述这3个算例的计算结果。

4.1 准静态模拟计算结果

第一个算例为准静态键型PD模型算例，计算对象为二维各向同性轴向受拉板(图8)，无预制裂纹，板长和板宽均为50 mm，板厚设为单位1，密度ρ=8 000 kg/m3，弹性模量E=192 GPa，泊松比v=1/3。离散化参数：x方向和y方向的物质点离散间距均为0.5 mm，近场范围为3倍的物质点间距。边界条件：在板的上下两端虚拟边界层处施加大小分别为2 000 mm/s和 -2 000 mm/s的速度边界条件，最大拉伸位移值d=0.2 mm，虚拟边界层厚度为3倍的物质点间距。总的物质点数为10 600，最大迭代步数设为5 000，时间步长dt=1.0×10-7s/step。为了模拟板在上下两端受拉情况下的弹性变形过程，不考虑物质点的损伤情况，将临界伸长率s0设为1。

图8 二维轴向受拉板

为了验证程序计算结果的正确性，分别取y坐标为25.25 mm处的一排点x方向和y方向的位移值，和理论解作对比，对比结果如图9所示。

图9 部分物质点位移值PD解和理论解的对比

从图9可以看出，除了靠近边界处的物质点位移值有微小的差异(误差<5%)，这是由PD数值模型自身所引起的，其余物质点位移值的PD解和理论解基本保持一致。表1给出了图8中标注的部分物质点在模型达到最终平衡状态时，x方向和y方向的位移值dispx，dispy，并将串行程序和并行程序的计算结果作对比，可以发现两者基本保持一致，从而验证了并行程序的正确性。

表1 串行程序和并行程序计算结果的对比

4.2 动态裂纹扩展模拟计算结果

第二个算例为二维各向同性矩形板裂纹分叉模拟算例。板长和板宽分别为100,40 mm，在板左侧中心处有一条长度lc=50 mm的预制裂纹。板的上下两侧有均匀拉伸荷载σ=12 MPa，如图10所示。材料参数：弹性模量E=65 GPa，泊松比v=1/3，密度ρ=2 235 kg/m3。离散化参数：x方向和y方向的物质点离散间距均为0.2 mm，近场范围为4倍的物质点间距，总物质点数为103 000。时间步长dt=5.0×10-8s/step，迭代步数设为1 000，临界伸长率s0=0.001 6。

图10 轴向拉伸下的单边预制裂纹二维矩形板

图11给出了模型在迭代步数为500和1 000时的裂纹扩展情况。可以看出，裂纹从预制裂纹的右端开始逐渐延伸，到一定阶段开始分叉。由于预制裂纹位置的对称性，裂纹扩展路径也呈现出完美的对称性(与已有文献[6]图8中的裂纹分叉结果基本一致)。

图11 算例二的计算结果

第三个算例为板中心裂纹动态扩展模拟算例，是在准静态模拟算例模型的基础上，在板的中心位置增添了1条预制中心裂纹，如图12所示，裂纹长度2a=10 mm。离散化参数：物质点离散间距为0.1 mm，近场范围为3.015倍的物质点间距，总物质点数为253 000。边界条件：在板的上下两端虚拟边界层处施加大小分别为50 000，-50 000 mm/s的速度边界条件。时间步长dt=1.3367×10-8s/step，临界伸长率s0=0.0447 2，迭代步数设为1 000，其他参数保持不变。

图12 轴向拉伸下含中心裂纹的二维板

表2给出了算例三在线程数设为1的情况下程序各模块的运行时间情况。

表2 算例三在线程数为1时程序各模块的运行时间

Insert center crack模块的作用是在板中心位置处添加预制裂纹，其余各模块的含义在前文已经介绍过，在此不再复述。由表2可以看出，在程序的整体运行过程中，相比其他模块，计算时间比较长的3个模块分别是Neighbor、Time integration和Compute kave模块。

图13给出了模型的初始状态，以及迭代步数为650和1 000时的裂纹扩展情况。从图13可以看出，裂纹从预制裂纹的两端开始逐渐扩展、分叉，形状犹如鱼的骨头(与文献[15]图11(b)、文献[16]图9.7中的裂纹分叉结果基本保持一致)，直到板完全被撕裂，整个动态裂纹扩展过程结束。

4.3 并行计算结果与分析

本文所采用的实验平台为Intel(R)Xeon(R)Core44 CPU E5-2699 v4，主频为2.20 GHz；操作系统为Windows10；开发工具为VS2012(支持OpenMP 2.0)。

上述3个算例程序整体运行时间与线程数之间的关系如表3所示，程序主要模块在不同线程数下的并行加速比情况如图14所示。

由表3、图14分析可知：随着线程数增加，程序运行时间呈下降趋势，且计算规模越大，粒子数越多，下降幅度越明显，得到的加速比越大，并行效率越高。在计算规模很小，物质点数不多的情况下，随着线程数增多，并行加速比并没有得到显著提高，而是趋于平缓，如图14(d)算例一所示。

图13 算例三的计算结果

表3 不同线程数下3个算例程序整体运行时间

这是因为运用OpenMP对程序进行并行化的过程中，线程组的派生和缩并，子线程的交互运行都需要花费额外的时间，这些都属于OpenMP自身的并行开销，而且随着线程数的增加，并行开销所占比例会不断增大，因此，当计算规模维持不变，程序运算负载不够大时，这些开销会影响到并行程序的加速执行，从而并行效率会降低。所以在问题规模固定的情况下，增加线程并不能持续提高计算效率，只有选择合适的线程数才能获得最快的计算速度[17]。因此，在实际并行化操作中，应根据具体的计算规模选择合适的线程数，避免因线程数量过多而速度更慢所导致的计算资源浪费和计算时间延长，或者计算规模应该随着计算线程数的增加而相应扩大。

从图14还可以看出，Neighbor模块的并行加速比情况比较理想，而其余几个模块的并行加速比没有显著提高。这是因为OpenMP主要针对无循环依赖的循环体进行并行化，如果程序中包含没有循环依赖的循环体，使用OpenMP能够大幅度地提高程序的并行性能。Neighbor模块中的循环体无循环依赖，而其余几个模块不仅有循环依赖的迭代循环，而且还存在子模块的调用，会增加额外的调用开销，因此，Neighbor模块的并行加速比比其他模块更为理想。

图14 3个算例程序中主要模块在不同线程数下的并行加速比曲线

5 结 语

在多核处理器平台上，基于OpenMP的近场动力学模拟相比串行平台，效率可以得到较好提升，随着线程数量的增加，加速比不断增大，计算时间大幅减少。这有利于充分利用计算机CPU资源，避免大量的CPU资源被闲置和浪费，也有利于近场动力学理论的进一步发展和应用。本文针对平面问题的键型PD模拟开展研究，今后需要进一步研究空间问题和态型PD模拟，在大型并行计算机上对并行算法和程序的可扩展性做进一步的测试与优化。