（计算数学专业论文）有限元若干问题的并行处理.pdf

摘要 有限元方法是解决当今各种大型工程分析问题的有效方法。随着规模的日益 扩大，有限元并行化也显得愈加重要了。已提出的有限元法并行化的集中途径： 硬件方面，在并行化时利用通用的硬件或设计专用的硬件，如各种向量机、有限 元机；算法方面，人们一直在研究并行算法，这些算法包括并行方程求解器及居 于区域分裂的并行算法、运算子分裂技术或单元接单元策略，以及在有限元分析 过程中的各个层次探索提高并行度的各种策略和技术。针对并行有限元的计算问 题，本文主要讨论了以下几个问题： l 、针对m i m d 并行计算系统，提出了一种特殊的待求解区域分裂方式 顺次分裂。这种特殊的方式使得并行计算系统在进行各子区域的网格划分时能充 分利用到现在的任何一种网格划分方法，不需要进行进程i 日j 通信。在顺次分裂基 础上又给出了特殊的结点编号方法顺次编号。按照这种编号方法形成的系统 方程具有特殊的形式。该形式保证了相邻子区域节点编号的统一，不用做任何数 掘传递，而且在求解系统方程时不合成总刚度矩阵。同时具有自适应的特点，可 以进行自适应有限元的计算。 2 、针对问题1 提出的由顺次分裂和顺次编号形成的系统方程，给出了求解 的多分裂并行算法；对经典的共轭梯度法实现并行化；利用多分裂法的分裂格式 构造了预处理共轭梯度法中的预条件矩阵，实现了预处理共轭梯度法的并行化。 同时还给出了算法的收敛性证明。 3 、利用问题2 中算法对两个偏微分边值问题进行求解，与串行算法进行比 较，对其并行性能进行了测试和相关分析。最后将算法应用于电磁场和弹性力学 的工程实际，说明了算法的有效性与实用性。 关键词并行有限元，网格划分，自适应，多分裂，共轭梯度法 a b s t r a e t f i n i t ee l e m e n tm e t h o di sa 1 1e 饿c i e n tm e t h o di ns o l v i n gt h el a r g ep r o j e c t n l e p a r a e lf i n i t ee i e m e n ti sb e c o m i n gm u c hm o r ei m p o r t a f l ta sar e s u l to fm ei a 唱e r s c a l e t h ee x i s t i n gw a y sa r e ：i n1 1 a r d w a f e ，u s i n gt h ec o m m o no rt l l es p e c i a lm a c h i n e m a tw a sd e s i g n e d ，f o re x 锄p l e ：e v e 叮轴n do fv e c t o rm a c h i n e sa n dm a c h i n e so i l i y f o rs o l v i n g6 n i t ee l e m e n t ；i na l g o r i t h m ：m a 芏l ys c h o l a r sk e e pd e v e l o p i n gt h e p a r a l l e l a l g o r i t t l i l l s t h e s ea l g o r i t h m si n c l u d i n gp a r a l l e le q u a t i o ns o l v i n gm a c h i n ea l 】dt l l e p a r a u e la l g o r i t l l i l lb a s e do nt h ea r e as p l i t t i n g ，t w o s t a g ec o m p u t i n gt e c h n o l o g yo r t h ee l e m e n tb ye l e m e n tm e t j l o d ，a n de v e r yk i n d so fm e t h o d sa n dt e c h n o i o g yi n e v e r ys t a g ed u r i n gt h ea n a l y s i so f t h e6 n i t ee l e m e n t t h i sp a p e ri sf o c u so nt h ef o l l o w e dq u e s t i o n s ： lp u tf o r w a r das p e c i a ls p l i t t i n g 、v a yo ft h es o l v i n ga r e ao nt h em i m p p a r a l l e l c o m p u t i r 唱s y s t e m _ s p l i m n gi no r d e l1 1 l i ss p e c i a lw a yc a nm a k eu s eo ft 1 1 ee v e r y e x i t i n gm e s hg e n e r a t i o nw h e ni tw a sm a d ei nt h ep a r a l l e lc o m p u t i n gs y s t e ma i l dw e d o n tm a k et h ec o 眦n u n i c a t i o nb e t w e e nt 1 1 ep r o c e s s e s as p e c i a in o d en u m b e f i n g m e t h o d m u m b 耐n gi no r d e ri sd e s i g n e db a s e do nt h es p i i t t i n gi no r d 托t h e s y s t e m i ce q u a t i o nf o r m e db yt h i sn u m b e r i n gm e t h o dh a sav e r ys p e c i a ls t r u c t u r e i t g u a r a n t e e sm es a m e n e s so ft h en o d en u m b e rb e t w e e nt h ea 锄a c e ma r e a w h e nt 1 1 e m e s hi sg e n e r a t i n gw ed on o tm a k et h ec o m m u n i c a t i o no ft h ed a 诅，w h e ns y s t e m i c e q u a t i o ni ss o l v i n gw en e e d n tt of o r mt h es t i f r n e s sm a t r i x t h i sm e t h o dc a i lm a l 【e t h ea d a p t i v ef i n i t ee l e m e n tc o m p u t i n g 2w eg i v et h ep a r a l l e lm u l t i s p l i t t i n gm e m o di ns o l v i n gt h es y s t e m i ce q u a l i o n f 0 册e di n1 ；m a k et h ec l a s s i cc o n j u g a c eg r a d i e mm e t h o db e c o m eap a r a l l e lm e t h o d ； f o mt h ep r e c o n d i t i o nm a t r i xo ft h ep r c c o n d i t i o nc o n j u g a t eg r a d i e n tm e m o db yu s i n g t h es p l i n i n gf o m l a to fm u l t i s p l m i n g ，a c h i e v et h ep a r a j l e l i z a t i o no ft l l ep r e c o n d i t i o n c o 玛u g a t eg r a d i e n tm e m o d 3w es o l v e dt w op a r t i a ld i f 五期e n t i a le q u a t i o n sb yt h em e t h o d sg i v e ni n2a n d c o m p a r e dw i mt h es e r i a lm e t l l o d s 1 1 1 e nw em a k et h et e s ta i l da 1 1 a l y s i st h ep a r a l l e l c a p d b i l i 印a “a s tw ea p p l i e dt h em e m o d si ne l e c t r o m a g n e t i cn e l da n de l a s t i c i t ya n d s h o w e dt h ee f n c i e n c ya n dp r a c t i c a b i l i t yo f t h em e t h o d s k e y w o r d sp a m l l e lf i n i t ee l e m e n t ，m e s hg e n e r a i i o na d 印t i v ef i n i t ee l e m e n t ， m u l t i s p l i n i n g ，c o 巧u g a t eg r a d i e n tm e t h o d i l 西北工业大学业 学位论文知识产权声明书 本人完全了解学校有关保护知识产权的规定，即：研究生在校攻读学位期间论文工作 的知识产权单位属于西北工业大学。学校有权鳔留弗向国家有关部门或机构送交论文的复 印件和电子版。本人允许论文被查阅和借阅。学校可以将本学位论文的全部或部分内容编 入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 同时本人保证，毕业后结合学位论文研究课题再撰写的文章一律注明作者单位为西北工业 大学。 保密论文待解密后适用本声明。 学位论文作者签名 j i i 明年月8 日 指导教师签名：墨丛 伊释j 月夕日 西北工业大学 学位论文原创性声明 秉承学校严谨的学风和优良的科学道德，本人郑重声明：所呈交的学位论文，是本 人在导师的指导下进行研究工作所取得的成果。尽我所知，除文中已经注明引用的内容 和致谢的地方外，本论文不包含任何其他个人或集体已经公开发表或撰写过的研究成 果，不包含本人或其他已申请学位或其他用途使用过的成果。对本文的研究做出重要贡 献的个人和集体，均已在文中以明确方式表明。 本人学位论文与资料若有不实，愿意承担一切相关的法律责任。 学位论文作者签名 年弓；j3 日 西北丁业大学颀卜论文 第一章绪论：有限元法的并行处理技术 1 1 引言 随着科学技术和工程技术领域的不断丌拓和发展，出现了各种各样大型和超 大型的计算的复杂结构。这些结构不但是具有很大自由度的组合结构，还包含非 线性本构关系、随机载荷和复杂的边界条件等多种因素。对这些结构的分析必须 借助于高阶的数值分析模型和大规模数学计算，以保证数值解的精度。而传统的 串行有限元分析方法受求解此类问题所需的计算时间长和内存储量大的限制，难 以实现这些结构的计算，使其应用受到束缚。因而需借助于已出现的各种各样的 并行机进行结构有限元并行处理。 目前己提出的有限元法并行化的集中途径包括：硬件方面，在并行化时利 用通用的硬件或设计专用的硬件，如各种向量机，有限元机；算法方面，由 于有限元结构分析中最耗时的计算阶段是重复求解大型线性代数方程系统。因 此，人们一直在研究并行算法，这些算法包括并行方程求解器及基于区域分裂的 并行算法、运算子分裂技术或单元接单元策略，以及在有限元分析过程中的各个 层次探索提高并行度的各种策略和技术【2 】。 1 2 有限元分析的并行处理技术 并行有限元分析理论与并行计算机的发展是密切相关的。早期的并行有限元 大部分都使用向量流水处理机、阵列处理机等专业并行机。自从美国国家宇航局 ( n a s a ) 的a k n o o r 于1 9 7 5 年发表第篇关于有限元并行计算的文章以来，有 限元并行计算的研究已经历了近3 0 年的发展历史，并取得了众多令人瞩目的研究 成果，包括有限元并行算法研究和并行程序设计等成果。程建钢【3 1 ，李晓梅4 】曾 对1 9 9 5 年以f ； 的并行有限元发展作了综述，主要讨论了并行有限元在向量机、阵 列处理机、共享存储多处理机等上的应用情况。 2 0 世纪9 0 年代，网络技术的快速增长对并行计算机的发展产生重大的影响， 例如，用专门网络紧密耦合多个计算结点而构成的大规模并行处理机系统大规模 并行机( m a s s i v e l yp a r a l l e lp r o c e s s o r ，m p p ) ，以及由多台高性能计算结点经计算 机网络连接组成的集群处理环境机群并行机( c i u s t e ro f w o r k s 诅t i o n ，c o w ) ，这些 典型的分布存储计算机系统逐渐成为高性能计算机领域的两大重要研究方向。尤 第一帝绪论：自限几的并 。+ 处世救术 其是c o w 具有投资风险小、结构灵活、可扩展性强、通用性好、异构能力强等 较多优点，已成为高性能计算领域的一个新的发展热点。 1 2 1 分布存储环境并行有限元的研究进展 在分布存储并行环境下，目前的有限元的并行策略主要有： 基于线性方程组并行数值求解器。对传统有限元最耗时的过程系统方 程组的求解进行并行化以加快整个分析过程的求解时间： 并行子结构法。根掘传统区域分裂算法的理论划分子域，在边界结点上建 立整体结构的平衡方程，降低问题的自由度： 并行e b e ( e l e m e n t b ye l e m e n t ) 法。在单元级上实施有限元整体运算的并行 化迭代算法，避免组集总刚矩阵； 其他方法，如f e t i ( f i n i t ee l e m e n tt e a r i n ga 1 1 di n t e r c o n n e c t i n g ) 法。 1 2 2 基于系统方程求解的并行有限元 并行有限元最直接和最仞的实现是利用数值并行求解器求解系统方程组，其 原理如图l ，1 所示，即单元分析和系统方程的形成是串行过程( 在一台处理机完 幽1 1 成) ，而系统方程组的求解是并行过程。采用这种方案的最简单做法是利用一些 p q 北t 业大学硕卜沦文 现成的线性方程组软件包。m a i l u e ll r d m e r o 【5 j 硬用s c a l a p a c k 软件包编制了 一个钢筋混凝土空问框架结构非线性分析的并行程序；李根国等【4 】将商业著名有 限元软件n a s n u n 在“神威i ”超级计算机系统进行了移植和并行计算功能的 二次开发。他们的丌发思路是将有限元的求解过程( 线性方稃组的求解或矩阵的 特征值问题) 进行并行化，并采用了p e t s c 软件包，前后的处理工作由n a s t r a n 软件完成，模块l 日j 通过文件进行数掘交换。 另一种做法是将刚度矩阵按行进行“分割”，各处理机存储与每块范围相应 的有关数据( 单元、结点、荷载等) ，这样单元分析和系统方程组装过程可并行进 行，不需任何数据通信( 当然是以少量的重复计算为代价) ，系统方程组的求解一 般采用迭代并行求解器( 如并行共轭梯度法，p c g ) 【。这种思想的解释如图1 2 。 原结构划分成若干子区域，划分面为“共享单元”( 图1 2 中的阴影部分) 。若每 个处理机保存各自子区域的相关信息( 包括“共享单元”的信息) ，则子区域内的 单元分析和总刚的集成是完全并行的，而只有“共享单元”涉及重复计算，计算 过程如图1 3 所示。 裘擎黪纛 霉 图1 2 第一审绪论：卡i 限儿的并 r 处理技术 进稃o 数据输入 建立有限元模刑 单元分析 集成系统方程 并行求解系统方程 更新单元状态 通信 例i 3 进榉p 数据输入 建立有限元模掣 单元分析 集成系统方稃 并行求解系统方稃 更新单元状态 当然，也可以采用基于传统的子结构划分方法( 即控制划分面是通过结点) 来实现如图3 所示的并行计算流程。不同的是，完整的系统方程需要通过通信、 叠加( 与边界结点相对应的部分) 后才能形成。实际上可直接进行求解，而不需组 装系统方程1 1 0 1 1 1 l 。显然，这种方案利用各“予结构”在单元分析和组装刚度方程 过程中的独立性来实施并行化。但是，整体有限元方程的阶数仍然是原结构的全 部自由度数，其求解时| 日j 仍占整个分析过程较大的比例。 总之，基于系统方程组并行求解的并行策略是一种较为简单、易于实现的方 法。由于系统方程组的系数矩阵是大型的、非常稀疏的，因此选择或开发并行求 解器时，必须要考虑到这一因素。此外，研究表明：对稀疏矩阵进行分解，若 仍采用与稠密矩阵相同的方法则并行效果很不理想，而需要采用有效重编弓和映 射算法j 能得到较好的效果；直接利用现成的并行软件包可实现并行有限元， 但是其效率有待充分评价。 1 2 3 系统方程组的并行求解 求解系统方程组的并行算法有两大类，一是并行直接算法，二是并行迭代算 法。从有限元分析的观点柬看，这两种算法都是在基本保持有限元法运算的框架 上对系统方程组实施并行计算。 并行直接算法是基于系数矩阵的直接分解，按分解的方式对它们进行分类， 4 两北厂业人学硕“仓艾 目前较多使用的算法有三分解算法、c h o j e s k y 分解算法以及稀疏线性方程系统 的向前和向后回代并行算法。并行直接算法的主要优点是算法稳定性好、精度高， 同时便于对现在广泛使用的串行算法程序j 兰行改造。其主要缺点是在方程的并行 求解过程中，由于多次调用系统功能束实现并行控制，浪费了不必要的进程等待 消息和传递消息叫问，导致并行算法的效率下降。 并行迭代算法是系统方程组并行求解的主要策略，并且自从上世纪6 0 年代早 期以来直在应用。迭代求解算法比直接求解算法要求的储存更小、精度更高。 从性能方面来讲，直接求解算法往往比迭代求解算法快并且在串行机上优先使 用。然而由于迭代求解算法固有的并行性，像预处理共轭梯度法( p c g ) 、最小残 余法、逐次超松弛法( s u c c e s s i v eo v e r r e i a x a t i o n ，s o r ) 、g a u s s s e i d e i 法等等，因 此此类算法在并行求解中重新引起注意。迭代算法可分为两大类：经典法( 静态 法) 和非静态法。经典法好理解且较老，非静念法是近年来丌发的，通常较难实 现和理解，但效率较高。如多色s 0 r 法、并行s 0 r 法、对称线性方程组的标准并 行p c g 算法等等。并行迭代算法的主要优点在于可避免和减少同步控制，有些方 法如e b e 方法，可以在不形成总刚度矩阵的情况下进行求解1 2 1 。它的主要缺点是 不便于和现有的些成熟的大型软件接口。 1 3 本文的主要工作 充分发挥高性能并行计算机的潜在性能的关键是：依掘不同类型并行计算机 的结构特点进行并行算法的设计和并行程序的实现。对于分布式存储大规模并行 计算系统束说，合理地分和数掘且尽量减少整体通讯是关键所在，要特别避免在 处理机之间的网络上的长路径上进行大舰模的数掘传递。 本论文的研究分为三个部分。 l 、针对m 【m d 并行计算系统，提出了一种特殊的待求解区域分裂方式 顺次分裂。这种特殊的方式使得并行计算系统在进行各子区域的网格划分时能充 分利用到现在的任何一种网格划分方法，不需要进行进程问通信。在顺次分裂基 础上又给出了特殊的结点编号方法顺次编号。按照这种编号方法形成的系统 方程具有特殊的形式。该形式保证了帽邻子区域结点编号的统一，不用做任何数 据传递，而且在求解系统方程时不合成总刚度矩阵。同时具有自适应的特点，可 以进行自适应有限元的计算。见第二章。 第一审绪论自限几的并行处刚技术 2 、针对问题1 提出的顺序分裂和顺序编号，给出了求解系统方程的多分裂 并行算法；对经典的共轭梯度法实现并行化；利用多分裂法的分裂格式构造了预 处理共轭梯度法中的预条件矩阵，实现了预处理共轭梯度法的并行化。同时还给 出了算法的收敛性证明。见第三章 3 、利用前面给出的算法对两个偏微分边值问题进行求解，与串行算法进行 比较，对其并行性能进行了测试和相关分析。最后将该方法应用于电磁场和弹性 力学的工程实际，说明了算法的实用性。见第四章 最后，对本文作出了工作总结与将柬工作的希望。 本论文是基于l i n u x 操作系统下的m p i 软件，在西北工业大学高性能计算 研究与发展中心的h pr x 2 6 0 0 并行计算系统上运算完成的。所有串行程序均是 用f o r t r a n 语言编写，并行程序是用在串行语言基础上扩展的并行语言 f o r t r a n + m p i 编写。 6 西北t 业人学坝l 论文 第二章有限元的网格并行生成 近1 0 年来，有限元网格生成方法研究有两个显著特点：( 1 ) 与其它研究领域 一样，经历了一个进化过程，一些方法的研究与应用出现停滞，而另外一些方法 在不断地深入、完善和发展，成为适应性强、。应用范围广泛的通用方法；( 2 ) 领 域和主题在不断扩展和深入，研究重点由二维平面问题转移到三维曲面和三维实 体问题，从三角形、四面体网格自动生成转移到四边形、六面体网格自动生成， 在并行网格生成、自适应网格生成、贴体坐标网格生成、各向异性网格生成等方 面亦取得许多重要进展。 2 1 有限元网格生成技术概论 有限元网格生成的通用方法主要有：( 1 ) 映射法，包括保角映射法( c o n f o n l l a l m 印p i n gm e t h o d ) 【4 16 】基于偏微分方程法( p d e - b a s e dm e t l l o d ) 【1 7 2 2 】以及代数插值 法( a l g e b r a i ci n t e r p o l a t i o nm e t h o d ) 【2 3 也6 】三大类；( 2 ) 基于栅格法，又可分为正则栅 格法( r e g i l l a r 鲥dm e t h o d ) 和有限四( 八) 叉树法( f i n i t eq u a d t r e e 、o c t r c em e t h o d s ) 两大类；( 3 ) d e l a u n a y 三角剖分方法( 简称d t ) 是目前最流行的通用的全自动网 格生成方法之一，大体上可将d t 算法分为三大类：分治算法、逐点插入法和三 角网生长法。 并行化计算环境对于大规模、超大规模科学计算以及高端工程应用是必需 的，而分布式计算环境可作为一种高端工程应用解决方案。现有网格生成并行化 或分布化算法在计算效率、内存管理、生成单元质量等方面还不够完善，还有许 多潜力可挖。另外，并行计算环境与分布式计算环境的控制软件日趋成熟，这为 算法的并行化、分布化丌发提供了更强有力的技术保障。先进的可伸缩的并行计 算机的出现为求解超大规模问题提供了可能性。有限元分析求解器的并行化已取 得了许多重要进展，但网格生成并行化研究却相对落后。当问题的规模达到百万 结点量级时，单机或串行机有限元网格生成就成了严重问题：一方面网格生成时 间过长，另一方面单机内存容量也不能满足要求。网格生成已成为大规模及超大 规模科学计算的瓶颈问题，因而网格生成算法的并行化是一个相当迫切的研究课 题。 衡量网格生成并行算法的主要指标是算法的可伸缩性、并行效率和网格质量 稳定性，其中最重要的是可伸缩性，因为它保证了当处理器增加时可在一个合理 第一审有限几的嘲格片行生成 时河内求解更大规模的问题。目前的网格生成并行算法都莎及到“分区( d o m a i n p a r t i t i o n i n g ) ”概念，所谓分区就是将目标域分解为多个子区域并提交给多个处 理器，好的分区应保证处理结点的负载平衡。目前的并行网格生成算法大体可以 分为三类【2 7 】：( 1 ) 子区域界面网格生成在i j i ；( 2 ) 子区域界面网格生成在后：( 3 ) 子区域网格生成与界面网格生成同步。 第一类算法是将目标域分割为多个子区域，并在网格生成之前预先生成界面 网格，这样在子区域的网格生成过程中就无需处理结点之间的通信。该类方法可 按予区域的分区形式进一步地细分为三个子类：初始粗网格分区，背景树分区， 表面网格直接分区。该类算法的优点是子区域网格生成不需要处理结点之丑j 的通 信，能够保障可伸缩性：缺点是由于在区域内部人工地加入了界面可能使网格质 量降低。 第二类算法是首先进行子区域网格生成，而子区域界面的网格生成曰待以后 处理。这类方法首先由s h o s t k o 和l o h n e r 【2 8 】提出，d ec o u g n y 和s h e p h a r d 【2 7 ，2 9 j 以及 张明敏和潘志庚等1 3 明提出的并行网格生成方法也属于这一类。该算法具有容错性 和自动负载平衡的能力。这类方法不需要处理结点之间的通信，能够保障可伸缩 性，然而再分区操作可能会降低算法的整体效率。 第三类算法与传统的做法不同，子区域网格生成与子区域之间界面的网格生 成同步进行，由c 嘶s o c h o i d e s 等【3 l 】提出的b o 、v y e r w a t s o n 算法的并行化实现是这 类算法的典型实例。这类算法的优点是在目标域内部未加入人工界面，所得到的 网格质量较好；缺点是需要处理结点之问的通信，不能够保障并行算法的可伸缩 性。此类算法有较大的发展潜力，通过对算法的精细设计甚至可达到比传统算法 更高的效率。 2 2 网格划分并行计算 结合第一章提到的基于系统方程求解的并行有限元，本论文提出一种特殊的 区域分裂方式和编号方法，使得划分网格和结点编号都可以并行计算，大大提高 网格划分的并行度。 2 2 1 子区域顺次分裂和结点顺次编号 针对以上算法的缺陷，本文提出如下的分区和结点编号方法，以矩形区域( 图 1 ) 为例，步骤如下： 西北丁业大学硕l 论文 ( 1 ) 将定解区域g 进行如图2 1 的划分，每个进程便区域必须两两相邻，使 其形成p 个子区域g 。、g 2 、g 。，每个区域作为一个不可分割的计算单位， 称为任务【l 】。每个处理器负责一个任务的计算。 g l g 2g 。 图2 1 域的划分 在划分时须注意以下两个原则： 子区域面积相等原则：即尽量使所划分的子区域的面积大致相等，这样有 利于每个处理器的负载平衡； 子区域两两相邻原则，即每个子区域只与前后两个子区域相邻，而与其管 子区域不相邻。 ( 2 j 在以上子区域划分的基础上，在每台处理器中利用2 1 节中提到的任何一 种网格划分方式进行并行网格划分，以四边形单元为例，可以得到如图2 2 的网 格划分。 l6 l l1 6 2 7 g 2 1 2 1 7 g g ， 3 8 1 38 4 9 4 图2 2 予犀域结点编号 ( 3 ) 在每台处理器内对结点进行结点编号。编号时须注意一个原则：每个子 区域与上一个子区域相邻的结点顺次编为这个子区域的最前的结点编号，而与下 一个子区域相邻的结点按顺次编为这个子区域的最后的结点编号。 9 第一帝青限儿的列 并并行竺成 按照以上的区域分裂方式和结点编号方法可以在每个处理器内形成特殊的 子区域刚度矩阵和子区域右端向量，其形式如下 眉”置篁k p 眉“= j 眉甜置封面必l ， ( 2 1 ) i 眉：眉掣酬，f m o 1 ) = l 。 ( 2 2 ) f f 其中上标f 表示第f 台处理器且，= 1 ，p 。而纠 、z 。是与上一个处理器的 重叠结点形成的子矩阵，叫，、。是与下一台处理器的重叠结点形成的子矩阵。 于是第一个处理器和最后一个处理器中的子区域刚度矩阵和子区域右端向量则 具有如下的形式： 北蹦甜嘲 ， 麟甜产) = 辫 由此形成的系统刚度矩阵则具有如下的形式： 置= 置磐石磐 d 置掣k 粤+ 眉丘g d 置! j 耳塞 d 丘置 ddd ddd 系统右端向量如下 1 0 l = dd j i g d 丘挣 d 眉掣+ 置0 d d 一o 1 ) + z ( 2 爿2 ) 一2 ) + ) f 擎一i 七f ， ”) do d0 dd dd 眉f + 眉0 置) 置身 眉2 ， ( 2 3 ) ( 2 4 ) 西北t 业人学帧l 论史 由于系统刚度矩阵和系统右端向量具有如此特殊的形式，所以在求解这样的 方程过程中就可以充分利用基于系统方程求解的并行有限元策略，使得子区域内 的单元分析和子区域刚度年阵的集成可以完全并行，而只需在求解方程组时进行 必要的数据交换即可。尽管这种算法存在共享单元必须进行重复计算的不足之 处，但是这些重复计算只占整体运算时间很小的一部分，完全可以忽略不计。 另外，当处理器台数确定时，该算法还可以进行自适应有限元的计算。其方 法可大致分为以下步： ( 1 ) 利用初始网格进行并行计算；若不满足精度要求，转( 2 ) ； ( 2 ) 在需要插入的位冒插入新结点，并对其进行总体编号且排在( 1 ) 中编号的最 后。对与新结点相邻的旧结点进行修改，将新结点的信息加入子区域刚度矩阵， 再将边界结点与新结点重新进行编号，以满足顺次编号原则，并置换矩阵记录新 编号与原有编号的关系。 ( 3 ) 该新生成的子区域刚度矩阵分别左乘、右乘置换矩阵。得到可以进行并行 计算的予区域刚度矩阵： 以两台处理器为例，求解区域为q ： ( 1 ) 划分初始粗网格如图2 3 所示。第l 台处理器中存储q 。的所有信息，其子 区域刚度矩阵如式2 5 所示。就此网格进行计算，若不满足精度要求，则转2 ： 置( 1 ) ： 图2 3 初始粗网格结点划分 t 开鹾 + t 芹 t 蠹+ 埒 。 簖 础 0 嗡+ 蠕k 譬 础+ 七箸埒 螺k 譬 ( 2 5 ) 第一章白限几的嗍格并行生成 ( 2 ) 假设需在图示位置加入一个新结点，由此形成新的单元如图2 4 所示。 图2 4 初始粗网格结点划分 这时由于加入的新结点5 将原来的单元化分为3 个新的单元，则只需要将原 租网格的单元的三个结点l 、2 、3 的信息进行修j 下即可，而不需要修改单元。 又因为引入了结点5 ，则只需要将原子区域刚度矩阵m 增加l 维，将第5 个结点 产生的元素存储在第5 行和第5 列即可，具体如式2 6 所示。 j r ( 1 ) ： 研+ 搿础 耘袅k 尝+ k 冀+ k 嚣 k 盆k 譬+ k 甓 。 瑶 k 基+ 碌k 翌+ 球 螺 螺+ k 雹 后岩+ 七譬+ 七譬 螺+ 蜷 而形成的置换矩阵如式2 7 所示。 p = olo o ol lo 0 o 础+ p 譬 蠕螺+ 螺 曙碍+ 七罢 碍 o q k 笔+ k 暑+ k 嚣 ( 2 6 ) ( 2 7 ) ( 3 ) 给新的子区域刚度矩阵茁o ) 分别左乘和右乘置换矩阵p ，得到下式。 主( 1 ) ：p 7 x ( 1 ) ，： 矧+ 臀 硝 o 蚓+ | j ；? 七；? 础 。 磋+ 七等+ 譬_ j 等 k 嚣嵫 + 硪 。 七；卜硪七筹 础+ p 等 女芸+ 礤 0 七要+ 后嚣+ 】 詈 碍+ 硪 瑶 罐+ k 嚣 瑶 t 器+ 磋 t 箸+ 蠼+ 磋 两北t 业人学坝i 论文 经过以上三个步骤，就达到了对局部求解区域进行加密的目的。 然而以上的算法有一个先决条件，就是处理器的台数必须保持不变。否则， 当处理器的台数发生改变时，每个处理器求解的区域就发生了变化，那么究竟在 那个处理器内进行加密就需要重新判断了。 第二市自限几系统方程的并行故值算法 第三章有限元系统方程的并行数值算法 3 1 有限元的子区域多分裂并行算法 由第二章的算法得到了以特殊的总刚度矩阵为系数矩阵的系统方程，下面利 用其系数矩阵的特殊性，为其设计了特殊的解法，为此先介绍求解线性方程组的 多分裂并行算法。 3 1 1 线性方程组的多分裂并行算法 并行多分裂迭代法是由0 l e a r y 和w h i t e 于1 9 8 5 年基于矩阵多分裂 ( m l l l t i s p l i t t i n g ) 提出的【3 3 】。其后许多国内外学者对此类方法进行了深入的研究， 使其方法得到了长足的发展。如b a i 【3 4 。6 1 ，b r u 等【3 7 q 9 1 ，f r o m m e r 等【4 0 4 ”。 求解大型线性方程组 a r = 6 ( 3 1 ) 其中，4 是非奇异矩阵。 一般的分裂迭代法可定义为： 腋o “) ： k “) + 6 ，i i ：o ，l ，2 ， 为引进并行计算机制，o l e a r y 和w h i t e i 叫提出了解线性方程组( 3 1 ) 的并行 多分裂迭代法，其基本思想是 令肘，、，和e ，都是矩阵，= l ，2 ，口，若满足 ( 1 ) 4 = 鸠一，圻1 存在；( 3 2 ) ( 2 ) e = j r ( j 为单位矩阵) ，日是非负对角矩阵则称三元组 f = l 帆，局y = 1 ，2 ，口为4 的一个多分裂。 把方程组( 3 1 ) 写成 x = m i j n t x + m 7 t bl = 1 ，2 ，一，a 用权矩阵e 组合这口个方程，得到 4 口口口 x = e | x = e | m - n | x + e 。m ? h = h x + g b ，= lk l，= l 西北工业大学坝r 论文 其中 口口 h = e | m _ n l ，g = e l m - 、 ，t l，s l 再对上面的方程组使用迭代 z o “) = j 勖“) + 劬，七= o ，】，2 ， ( 3 3 ) 此迭代格式称为多分裂迭代格式，日称为多分裂迭代法的迭代矩阵。于是产生 如下算法 多分裂迭代算法： 任取初始近似x ( 0 ) r ” 对j i = o ，1 ，2 ，重复做下面的( 1 ) 、( 2 ) 两步，直到收敛。 ( 1 ) 对，= 1 ，2 ，口，求解j ， 肘，j ，= ，工“) + 6 ( 2 ) 计算 x i ) _ 妻e j ， 下面我们来考察该算法的并行性开发。首先，我们看到在( 1 ) 中每个n 的计 算是相互独立的，因此，若有一台具有一个主机与口个处理器的并行机，则该算 法可这样并行实现：每台处理器处理一个局部迭代( 1 ) ，在每一步j 将局部迭代 值朋送到主机，主机对所( ，= l ，2 ，口) 进行加权平均而得到整体迭代值工( “) ， 即( 2 ) 的计算，然后将x ( 川) 送回各处理器以开始下一步的局部迭代。注意到如果 e 的某个对角元素为零，则乃的相应分量无需计算，从而大大节省工作量。由 此可见，局也起到调配各处理器上工作量的作用。因此也称局为权矩阵。我们 应尽量选择e ，使得各处理器问大致做到负载平衡，从而减少同步等待的丌销。 可见，该算法具有自然的并行性。 3 1 2 子区域多分裂并行算法( f m ) 以p o s s i o n 方程 第二市自限几系统方程的j 4 一丁数值算法 一窘一窘厂b ) e g ( 3 4 ) 1 。笺矿z ：誊。 。4 其中置、分别如( 2 3 ) 、( 2 4 ) 式所示。每台处理器内的子区域刚度矩阵和子区域 ，一豳 b s ， 对于第f 台处理器，由于其所储存的系数矩阵为眉。利用多分裂算法进行 k 卜 ， ， ( ，) ：) 一置： 而 1 6 防硎 0 ) ：l 置g l 础， 删置 眉封置粤l ， 磐置纠 d + 一1 ) 磷) 一砖”d d ddd dd 一1 ) 列一贰“ 一) _ 帜( 。，三，j ，圭加) d ，( 3 8 ) 西北t 业大学硕论文 即e o j 与总刚度矩阵同阶，只在置“ 对应的位旨有非零元，且在对应位置上子区 域与其它区域重叠的位置均取三，非重叠位置取l 。 以三台处理器为例，系统刚度矩阵如下： 置= 置髫置磐 砖p 鬈磐+ 置p d 置! ： d j j ：1 dd ddd 置g 鬈g o 置! ； 眉 o 置五+ 置g 置g d 置茁 各处理器中的肘( “、 r o ) 和e o ) 分别如下： 肘( 1 ) ： - ) ： e ( 1 ) ： 肼( 2 ) ： 胃粤 置身 ddd l i 置劣置9 ddd i d口 量! ；) ddl l dod 五+ 耳搿 d i 【d ddd 耳塞j ddddd d0 1 ) 砖2 一耐2 ) 一量翟一置冀d o 一磷。一磷 o o x 袋一k 袋。一x 毽 oook 娥 。 抄 置磐 口 d 置甲 d 置 ： 口 五j ； dd d d d d 置g 鹂) 置碧 d dd 置g d 置g o p 嚣g d d 五曼 1 7 第一章柏限几系统方程的并行数值算法 ( 2 ) ： e ( 2 ) ： 肘( 3 ) ： ( 3 ) ： e ( 3 ) ： d 抄 世 三舻 置磐 dd d 尉，+ 置掣 d dd 眉塞 dd口 ddd d 一置磐 d 一置磐 d 一置p d 一石l ； d d 一置笋一眉； dod d d 抄 d d dl odl doi 、 置舻量2 j 眉掣置塞j dd 一置 d 面! ； d 0 1 ) ，眉 ；一；d dd 由于权矩阵e 的设置，在进行并行求解系统方程组时，每台处理器只需要 与相邻的处理器进行数据交换即可。这样就避免了大规模数据传送的问题，减少 了通信的时间，提高了并行效率。由于对方程m ( ) j ，( ) = ( ) “) + 厂) 的求解比 较复杂，因此这罩采用定常二级多分裂迭代，即在每台处理器内进行方程组求 解时，令三( ，) 为膨o ) 的严格下三角矩阵，d ( ) 为m ( 的对角矩阵，【厂o ) 为肘( ) 的严 格上三角矩阵，使用与j a c o b i 算法相应的g a u s s s e i d e l 迭代算法迭代c ( c 为一常 正整数) 次求得j ，o ) 的近似值。于是形成有限元多分裂并行算法( f m ) 如下： 1 8 o d d 畦。 一 i 置一 d o d 蚶础 卜一 一 九铲 d d d d d m 驺 f一 砖卿d d d 一卜 一 m ” d 置d d d 两北t 业人学硕l 论文 ( 1 ) 在每台处理器中，同时进行如下步骤： 对每个单元，计算单元刚度矩阵置。和右端荷载向量。； 在每个处理器内合成子区域刚度矩阵置( ) 和子区域右端载荷向量，( ，并 施加边界条件； ( 2 ) 在每台处理器中重复做下面的( 1 ) 、( 2 ) 两步，直到收敛。 对扛1 ，2 ，求解j ，o j ： 6 0 ) = e ( ( ) 洲+ e ) ) ； j ，( m ) = ( i ) ) _ 1 ( - ) j ，( 枷u ( ，) ，) + 6 ( ) ，= o l ，2 ，c ，c 为一常正 整数； 计算 忙“) _ 艺e j 。 3 1 3 并行步骤 求解有限元的步骤如下，其中第一上标f 表示第f 台处理器，第二上标七表示 第忌次外迭代： ( 1 ) 在每台处理器中同时计算单元刚度矩阵眉8 和右端荷载向量，8 ； ( 2 ) 在每个处理器内同时合成子区域刚度矩阵置和子区域右端载荷向量 ( “，第f ( f = 1 ，p 一1 ) 台处理器取得第f + l 台处理器的 ( ) ，第f + 1 台处 理器取得第f 台处理器的) ；第f 台处理器计算j = j + z ) ，第p 台处理 器计算 ( ，) = 一川) + z ( p ) ；施加边界条件 ( 3 ) 采用多分裂法求解方程组： 取初始向量“( ； ， 每台处理器令“) = ，( f ) ，第f ( f = 2 ，p 1 ) 台处理器分别计算 p f ) = k 磬，置篓，置等- 卜1 和v ) = ，置，眉b ，- “+ ”，并将v f 传递至第f + 1 台处理器，v 传递至f l 台处理器；( 注：第1 台处理器只计算和传递v f l ) 第 1 9 第二审存限几系统程的j 仃数值算法 p 台处理器只计算和传递v 。) 第f ( f = 2 ，p 1 ) 台处理器接受第f 一1 台处理器传递的p p ) 和第f + l 台处理器传递的v p ) ，并计算z f ) = ( 。) 一v p 和= 爿“一v p ) 。( 第1 台处 理器只接受v 9 和计算”，第p 台处理器只接受p p 。和计笋五如) 第f ( f ：2 ，p 一1 ) 台处理器计算v f ) = 姐一1 ) 置”，d k 咄。”和 v l ) = ( d ，缸一1 ) 置磬- “”，并计算z = “一v 7 和爿2 = 一“一p p 。( 注：第 1 台处理器只计算p 9 和”，第p 台