第二十三届全国高层建筑结构学术会议论文2014年 并行计算在弹塑性动力时程分析中的实现和效率 刘慧鹏，李修宇1 （1.广州数力工程顾间有限公司，广州510000) 提要：本文闸述了并行计算在弹塑性地震动力时程反应分析中实现以及优化方法，采用CPUGPU的异构平台，重点闸明并 行计算算法、数据结构、显存利用以及优化方法等方面对SAUSAGE软件计算效率的作用，表明对弹塑性计算效率提高非常 显著.
关键词：并行计算，优化，弹塑性动力时程分析 1前言 弹塑性地震动力时程反应分析对建筑结构分析和设计都是非常重要的，它一定程度上能够反应地震作 用下的真实运动情况”.
但是，弹塑性地震动力时程分析就不可避免的遇到一个间题，就是对运动方程的 求解： M{] C[6] K{6]: = Pt (1) 式中[)、}、{6为某时刻加速度、速度、位移.
M、C、K为质量矩阵、阻尼矩阵、刚度矩阵，下同.
这个大型矩阵方程的求解是个高消耗、高难度的过程.
运动方程的求解大体分为隐式方法和显式方法.
但是这两种方法都不可避免的需要高耗时、耗资源的求解过程.
隐式求解方法是指用当前步结果和下一步 未知结果反复迭代下一步结果，必须通过选代得到.
显式求解方法指用上一步的结果和当前步的结果计算 下一步的计算结果.
国内从上世纪60年代就有很多人开始对运动方程的求解进行研究，并出现了很多对 运动方程求解的方法，但是总体运用到工程实际中的应用研究却非常缓慢.
虽然国内外近几年随着计算机 的发展出现了EPDA、SAP、PERFORM-3D等能够进行地震动力时程反应分析的软件，但是在精细化有限单元 和数值计算方面没能够充分发挥计算机硬件的优势.
并行计算是相对于串行计算而言的，指用多个处理器并发的执行计算，是提高计算机系统计算速度和 处理能力的一种有效手段.
并行计算的发展在近些年发展快速，尤其随着硬件在摩尔定律促进下快速成长 的同时，并行的概念也逐步深入各个行业当中，目前在图形、大数据、航天、大气气象等领域已经得到充 分的应用，但在国内工程结构领域里却未能得到充分应用.
在广州建研数力建筑科技有限公司的努力下，采用显式积分算法结合CPUGPU异构平台并行计算很大 程度的解决了运动方程求解过程的困难，使弹望性地震动力时程分析方法能够完美的应用到实际工程当中.
2SAUSAGE中显式积分算法的解耦 运动方程是个大型矩阵的求解过程，其中M、C、K均是N*N的大型方阵，跟节点自由度有关，考虑到 充分利用并行计算需要对这个大型方程组进行解耦成单个自由度下的方程.
因此需要对质量矩阵、阻尼矩 阵和刚度矩阵进行解耦.
首先将运动方程按照中心差分进行代换，即得到： (Z) 式中{8}ar、{8]、{6-ar分别表示下一时刻、当前时刻和上一时刻的位移，△t为时间步长，P为当前时 作者美介：刘草鹏（1983-)，男，张士.
国家一级注册结构工程师
第二十三届全国高层建筑结构学术会议论文2014年 间荷载向量，质量矩阵通常可采用集中质量对角矩阵，可以将单元的质量通过形函数集中到节点上，形成 对角矩阵.
阻尼矩阵一般可采用其他方法来计算，在SAUSAGE中可以选瑞雷阻尼和振型阻尼.
瑞雷阻尼即 把C分解为aMBK，当取β=0时，阻尼矩阵C也跟质量矩阵一样变成对角矩阵，求解计算量大大减少.
而a与结构周期成反比，有α=4punt/T，其中u根据规范取相应的阻尼比，T为结构周期，通常取结构第一周 期值.
振型阻尼则通过振型空间代换求得广义阻尼矩阵，即 [μM 1 C=0c-= (3) μMJ 在通过构件矩阵求逆可得到阻尼矩阵.
面K乘以位移则用等效不平衡内力来确定，先在单元内求出单元内 力再组装到节点上形成为节点内力向量.
这样既可实现方程的解耦，避免了刚度矩阵组装，矩阵求逆等耗费大量时间和资源的工作，而且解耦 后的方程可以充分应用到并行计算当中，可以做到每个线程求解一个方程组，每个线程进行一次等效不平 衡内力求解.
3并行计算的实现 并行计算的基本思想是对同一个问题分解成多个子问题，然后交给多个线程分别求解子问题.
在上一 节中把方程解耦成多个子方程，然后根据并行计算思想把每个子方程交由多个线程来同时计算，如下图所 示： 任务 任务 手任务1 承任务2 子任各2于任务1子征务3 承任务3 图1串行和并行 在并行计算的实现中最重要的要点在于数据结构的合理性和显存的合并访间要求.
sausage 中针对并 行计算的特点建立了自由扩展的数据结构.
在有限元计算中，首先要建立节点、单元、荷载数据结构，由 于荷载已经全部导算到节点上了，可以将质量和荷载都依附到节点信息内，并形成合并访问数据格式，并 且将节点的编号都按照顺序排列，即用数组下标做为节点编号.
如下： struct Node double* pCoord; //节点坐标 double* pMass; 节点质量 double* pLoad; //节点荷载 } 单元上需要依附截面序号、材料信息以及单元中的钢筋属性，需要定义钢筋的纤维束或纤维层的面积 和位置.
其中二维单元采用分层壳模型，采用统一划分层数，即在墙板单元信息中附带单元厚度即可.
struct Elem int* pNodeID; //单元节点 int* pSectID; //单元截面序号 int* pMatID; //材料序号 double*pFiberlnfo;//单元钢筋纤维信息
第二十三届全国高层建筑结构学术会议论文 2014年 } SAUSAGE中梁单元采用的是纤维束模型，因此可在CPU端提前划分好截面的纤维束，并按照纤维束的 面积和局部坐标存储形成合并访问数组.
在显卡中有多个core，每个core上有多个pe，每个pe计算一个方程过程，但是在任何硬件上core 和pe的数量是有限的，比如目前高端的nvidia显卡GTX780只有192个核心1024个pe，也就是每次进入 显卡计算的只有1024个子过程，如果每个pe上计算时间几乎一致，那么就这1024个pe都得到充分利用， 否则会出现木桶理论的短板，则这个批次中的1024个pe的计算时间是最慢的那个pe计算时间.在sausage 中影响每个pe计算时间的主要是纤维的个数和本构的一致性，也就要求每次进入核心计算的是同一个单 元类型和同样的截面类型.
因此根据建筑结构的特点区分主梁、次梁、杆元、边缘构件、墙和板构件形式 多个批次进入计算，每种构件类型具有特有的纤维束或纤维层数量，这样基本保证每次进入计算的多个线 程用时几乎一样的，达到负载平衡.
如下图2所示 主梁 36根纤维束 次操 12根纤维束 结 杆 1根纤维束 构 边缘构件 1根纤维束 培 6层纤维层 板 2层纤维层 图2结构分构件批次对应表 4并行计算的优化 计算效率的关键在于并行计算的效率，并行计算并非简单的实现即可，在定义良好高效的数据结构和 组织流程外，还需要对并行计算部分进行优化.
良好高效的数据结构可以减少通信，多批次分批进入计 算的区分可以使得有良好的负载平衡，这两点是并行计算优化的重点，在上一节中详细描述了.
但是在并 行计算细部而言也有需要注意的，而且对并行计算效率有很大的影响.
首先是对寄存器的使用数量，它是用来暂存指令、数据和地址的，能够存储的数量是恒定的，也是所 有存储设备里速度最快的部分，但当前活动的变量超过其存储的限额时就不得不转换一部分到内存中，即 出现寄存器溢出，由于内存的速度要比寄存器速度慢，显然当出现寄存器溢出时，并行计算的速度受到影 响，而且异出部分是哪个变量是随机的，如果溢出的是使用频繁的变量，必然导致并行计算效率大打折扣.
因此，SAUSAGE中对寄存器使用进行优化，尽量减少寄存器的溢出.
在并行计算中尽量少使用逻辑判断.
因为SAUSAGE的并行计算是在GPU中进行的，GPU对逻辑判断 非常不擅长，会影响计算速度.
而且当出现逻辑分支时，由于分支处理可能出现较大的差别而导致负载不 平衡，导致一些线程悬挂.
如下所示： else X=... x=iSTrae*(.)(1-srue)*() end if X... 图3逻辑分支转换 对加减乘除的使用也是值得优化的.
通常除法的效率是最低的，尽可能的少用除法，改用乘以倒数.
此外，优化四合运算的次数，将一些仅算一次的常量计算部分可先在CPU端计算完并存储起来，逾免反复
第二十三届全国高层建筑结构学术会议论文 2014年 计算.
5工程效率对比 SAUSAGE软件目前已经在全国大面积使用或试用了，其正确性和可靠性得到了相应的验证，这里只 讨论计算实际工程的计算效率，并与国际某通用有限元软件相对比.
如下图4所示工程， 图4模型三维简图 共8万多单元节点，7.6万多壳元，50万自由度，计算一条30秒长的地震波，计算步长为1.02e-004，采 组成的服务器上计算用时4个小时零49分.
如莊SAUSAGE中为了考虑阻尼更全面点，而采用振型阻尼， 随着采用振型数的增加而线性增加，如下图5所示： 300 250 200 150 100 50 时间 0 0 10 20 30 40 50 60 图5振型-计算时间图 上图横轴为振型数，竖轴为计算时间，当振型数为0时表示采用瑞雷阻尼.
从上图可以看出，如果采用振 型阻尼则当采用振型数为50个以上才与国际某大型通用有限元软件采用瑞雷阻尼计算的时间接近.
下表1 置和上面的一样.
表1计算效率对比表 序 工程名称 结构类型 层数 基本周期计算时间 某软件时 号 (h) 间（h) 1 顺德保利商务中心 框架核心筒 47 4.48 2.4 18 2 东莞长安万科中心 框架核心筒 60 6.5 2.8 19 3 成都世茂猛追湾（一期）8# 剪力墙 48 3.01 2 12
第二十三届全国高层建筑结构学术会议论文 2014年 塔楼 1 4 郑州华润中心二期5号楼 剪力墙 56 3.34 3.2 25 5 青岛华润中心悦府一期 剪力墙 63 4.43 3.6 24 越秀星汇云锦商业中心A 6 区A栋 框支剪力墙 50 3.72 3.2 19 7 华润惠州小径湾酒店 框架剪力墙 11 1.08 2 12 8 成都西部金融中心 框架核心筒 57 5.58 4.4 29 9 青岛华润中心悦府二期 剪力墙 67 4.3 3.6 23 10 华润深圳湾住宅 框支剪力墙 46 3.78 2.4 15 11 成都世茂猛追湾（一期）9# 塔楼 剪力墙 55 3.01 2 12 12 杭州华润MT楼 框架核心筒 60 2 15 13 成都顺江路333号B塔 框架剪力墙 67 6.14 4 29 14 山和华商贸广场 框架核心筒 40 3.48 3.6 24 15 渤海银行业务综合楼 框架核心筒 52 4.89 4 26 16 天津富力城 框架核心筒 92 7.31 5.2 32 17 郑州华润中心二期6号楼 剪力墙 55 3.78 3.2 21 18 天津现代城酒 框架核心筒 48 4.98 3.6 22 6结论和展望 SAUSAGE中采用并行计算来进行非线性弹塑性地震反应分析，在解耦运动方程后，通过数据结构的设 计、合并访问、负载平衡等方式大幅提高了计算效率，同时对并行计算进行局部的优化，使得计算弹塑性 分析实现在个人PC机运行的可能.
经过实际工程验证计算可行，而且比当前一些国际通用有限元软件的 显式计算模块计算速度提高了几倍.
随着硬件不断的发展，SAUSAGE将进一步的优化并行计算，提高计算 效率.
参考文献 [1]汪大绥，李志山，李承铭.复杂结构弹塑性时程分析在ABAQUS中的实现[J].第十九届全国高层建筑结构学术会议论文集， 2006. [2]常磊，李志山.SAUSAGE 在超高层建筑结构罕遇地震弹塑性时程分析中的应用.建筑结构[J]，2012.11 42(S2）：12-17 [3]黄忠海，廖云，王远利，某复杂超高层结构罕遇地震弹塑性时程分析[].建筑结构，2011，41(3）：40-44 [4]李志山，容柏生，高层建筑结构在罕遇地震影响下的弹塑性时程分析研究[J].建筑结构，2006，36（S1）：142-149 [5]OpenCL编程指南[S].机械工业出版社2012
