（计算数学专业论文）多核计算机集群上的有限体积元自适应并行计算.pdf

中文摘要 n a v i e r - s t o k e s 方程描述的是非常普遍的流体运动现象，因而具有十分广泛 的实际应用。这种方程的解通常具有异常大量的高频细节，因而其快速 数值求解常常需要用到自适应并行计算。有限体积元方法( f v e m ) ，又 称广义差分方法( g d m ) ，具备很多适合计算此类方程的性质。本论文 因此研究运一方法的自适应并行计算，特别是其自适应网格在并行计算中 的实现方法。 本论文在研究过程中，发现近年来出现的多核处理器对并行计算的主要 形式也就是集群计算有极大的负面影响。尽管多核处理器是单核 处理器的盛然继任者，尽管它是未来处理器的主流发展方向，尽管它具 备更高的本地计算能力，多核处理器非但不会顺其自然地提高计算机集群 的性能，而且还很可能使很多原有的并行计算程序变慢。美国加州大学 伯充菜分校的系统研究，以及美国国家实验室最近的模拟实验，都印证了 这一观点。宄其原因，其奂面作用主要源自于一个事实：多核处理器让 多个计算单元共同使用同一计算节点上的网络接口控制器，从硬件上大幅 度地增大了网络接口控制嚣的工作负荷。原有的并行程序没有清醒地认识 刭这一事实，自然就无法避短扬长，真正开发出多核处理器的计算潜力 了。 本论文因此而研究：在多核计算机集群上的有限体积元自适应并行计算 中，如何去解决这个网络负荷的问题。由于多核处理器是新生事物， 这 方面的研究目前也只是刚刚开始。本论文在资料搜集中只发现两类处于探 索阶段的解泱方案：( 1 ) 对计算机集群的硬件、操作系统和计算程序进 行彻底的系统的改造二这一方案尚未实现；( 2 ) 重新设计程序中的并 行算法，让其算法显式地利用多核集群计算中的多层内存布局这一方 案实现了一些基础算法的重构，包括稠密矩阵乘法和快速傅利叶变换的并 行计算。第一类方案的可行性有待考证，第二类方案在本论文的计算问 题中遇到了困难有限体积元的网格调整过程丛须使用保持网格质量的 自适应算法，其算法很难显式地为多层内存布局重构。因此，本论文另辟 新的途径，转而重构并行计算下面的基本通信模块( 数据交换、并行求和 i 以及并行求前缀和) ，从而即使保持原有的并行算法不变，都可以避开 多核处理器对集群计算的负面作用，证计算性能随着核数的增加而真正得 到提高。 本论文的论述主要分为四个部分：( 1 ) 解释有限体积元网格自适应 需要保持网格质量的原因和细节；( 2 ) 建立数据交换在多核集群计算中 的数学模型并且重构数据交换；( 3 ) 重构并行求和以及并行求前缀和； ( 4 ) 分析有限体积元的网格自适应并行算法在数据交换重构之前和之后的 性能变化。 本论文还通过实验证实了基本通信模块重构的可行性。实验数据还显 示，基本通信模块重构对当前的每节点四核的集群计算效果显著。因为以 后每过两三年多核处理器的核数都会翻一翻，本论文的方法值得继续观察 和进一步的研究。 英文摘要 n a v i e r - s t o k e se q u a t i o n sg o v e r nt h ew i d e s p r e a dp h e n o m e n ao ff l u i dm o t i o n t h e ya r et h e r e f o r ei m p o r t a n ti np r a c t i c e t h e i rs o l u t i o n su s u a l l yc o n t a i ne x t r a o r d i n a r i l yl a r g ea m o u n to fh i g hf r e q u e n c yd e t a i l s ，t h u st h e r ei sn e c e s s i t yt o r n nt h e i rn u m e r i c a lc o m p u t a t i o no na d a p t i v em e s h e sj np a r a l l e l f i n i t ev o l r i m ee l e m e n tm e t h o d ( f v e m ) ，a l s ok n o w na sg e n e r a l i z e dd i f f e r e n c em e t h o d ( g d m ) ，h a sm a n yf e a t u r e ss u i t a b l et ot h i sk i n do fc o m p u t a t i o n t h i st h e s i s t h e r e f o r es t u d i e st h ea d a p t i v ep a r a l l e lc o m p u t a t i o no ft h i sn u m e r i c a lm e t h o d ， e s p e c i a l l yt h ei m p l e m e n t a t i o no ft h ep a r a l l e la d a p t i v em e s h i n g i tw a sd i s c o v e r e dd u r i n gt h er e s e a r c hp r o c e s st h a tt h em a n y c o r er e v o l u t i o no fc o m p u t e rp r o c e s s o r si n c u r ss e r i o u sn e g a t i v ee f f e c t so nt h em a i n p a r a l l e lc o m p u t i n gf o r m c l u s t e rc o m p u t i n g t h o u g ht h et r a n s i t i o nf r o m m o n o - c o r et o m a n y c o r ei si n e v i t a b l e ，t h o u 曲m a n y - c o r ei st h em a i n s t r e a m o f f u t u r ep r o c e s s o rd e v e l o p m e n t t h o u g hm a n y c o r ep r o c e s s o r sp r o v i d ew i t h h i g h e rl o c a lc o m p u t i n gp o w e r st h a nt h e i rm o n o - c o r ea n c e s t o r s m a n y - c o r e p r o c e s s o r sc a n n o ta u t o m a t i c a l l yp r o v i d ew i t hh i g h e rc l u s t e rc o m p u t i n gp o w e r s ，a n dt h e ye v e nm a k em a n ye x i s t i n gp a r a l l e lc o m p u t i n gp r o g r a m ss l o w e r b o t hu cb e r k e l e ya n dt h en a t i o n a jl a b so fu n i t e ds t a t e sh a v ep o i n t e do u t t h i sf a c ti nt h e i rr e p o r t s t h et r o u b l e so fm a n y c o r ec l u s t e rc o m p u t i n gm a i n l y c o m ef r o mt h ec o n t e n t i o n st ot h en e t w o r ki n t e r f a c ec o n t r o l l e r ( n i c la m o n g m u l t i p l ep r o c e s s o rc o r e s t h ee x i s t i n gp r o g r a m sa t e n o ta w a r eo ft h i sf a c t ， t h u st h e yc a n n o tf u l l yu n l e a s ht h ec o m p u t i n gp o w e ro fm a n y - c o r ec l u s t e r s f o rt h e s er e a s o n s ，t h i st h e s i ss t u d i e sh o wt os o l v et h en e t w o r kc o n t e n t i o n p r o b l e mf o rf v e ma d a p t i v ep a r a l l e lc o m p u t i n go nm a n y c o r ec l u s t e r s s i n c e i t i ss t i l lt h ee a r l ys t a g eo fd e v e l o p m e n to fm a n y - c o r ec l u s t e r s t h ep u b l i s h e d s o l u t i o n st ot h et r o u b l e so fm a n y c o r ec l u s t e r sa r er a r e o n es o l u t i o ns u g - g e s t sam u l t p l a y e rr e d e s i g no fh a r d w a r ea n ds o f t w a r e t h eo t h e rs u g g e s t sr e - d e s i g n i n gp a r a l l e la l g o r i t h m sa n dr e w r i t i n gt h ec o d eo fu s e rp r o g r a m s t h e s e t w os o l u t i o n sa r eb o t ht o oh a r dt oa p p l yt op r a c t i c a lu s e e s p e c i a l l yt ot h e i i i a d a p t i v ep a r a l l e lc o m p u t i n go ff v e m w h i c hr e q u i r e sf o rt h eu s eo ft r e ea l g o - r i t h m st oa c h i e v ef v e m si n h e r e n tr e q u i r e m e n tt om e s hq u a l i t y t h i st h e s i s p r e s e n t sa n o t h e rs o l u t i o nt h a to n l yr e q u i r e st h ee x i s t i n gp r o g r a m st or e - c o m p i l e da n dl i n k e di nt h en e wp a r a l l e lp r o g r a m m i n ge n v i r o n m e n t ，i nw h i c h t h eu n d e r l y i n gb a s i cc o m m u n i c a t i o nm o d u l e sa ler e f o r m u l a t e df o rm a n y - c o r e c l u s t e r s t h em a i nb o d yo f t h i st h e s i si n c l u d e sf o u rp a r t s ：( 1 ) t h ee x p l a n a t i o n t of v e m si n h e r e n tm e s hq u a l i t yr e q u i r e m e n t ，( 2 ) t h ep e r f o r m a n c em o d e l o fd a t ae x c h a n g eo nm a n y - c o r ec l u s t e r sa n dt h er e f o r m u l a t i o no fd a t ae x - c h a n g e ，( 3 ) t h er e f o r m u l a t i o no fp a r a l l e lb u ma n dp r e f i x - s u m ，a n d ( 4 ) t h e t r e ea l g o r i t h m si nf v e ma d a p t i v ep a r a l l e lc o m p u t i n ga n dt h ec o r r e s p o n d i n ge x p e r i m e n t s t h ee x p e r i m e n t sp r e s e n t e di nt h i st h e s i ss h o wt h ef e a s i b i l i t yo ft h es o h t i o n t h e ya l s os h o wt h a tt h ep e r f o r m a n c ei m p r o v e m e n to ft h es o l u t i o ni s s i g n i f i c a n te v e no nc l u s t e r so fo n l yf o u rc o r e sp e rn o d e s i n c et h ec o r ec o u n t p e rn o d ew i l lb ed o u b l e de v e r y2t o3y e a r s ，t h es o l u t i o no ft h i st h e s i sa x e w o r t h yf o rf u r t h e rs t u d ya n do b s e r v a t i o n i v 原创性声明 学位论文作者签名： 日期：2 研年5 月 簇吱文争 2 日 ，容作文人 下内的在本 导的过已由 指用写均果 的引撰，结 师明或体律 导注表集法在经发和的 人已经人明 本中已个声 是文体的本，除集献到 文。或贡识 论果人要意位成个重全学的他出完 的得其作人 交取何究本 呈所任研。 所作含的明：工包文标 明究不本式 声研文对方 重行论。确 郑进本果明。 人立，成以担本独外口们中承 学位论文使用授权声明 学位论文作者签名： 日期：力孵f , e l 导师签名： 日期：之的7 年石月 许必i - 2 日l ， ：论量将缩 即交少权、 ，送的有印 定构的，复规机目阅用 的定利杏一采 文指赢被以 论其非室可 位或于料， 学门用资索 用部文系检使管论院行、主位、进留家学馆库 保国将书据。 关向权图数、叉有并有校关论 学文，学有位 大论版入入学山位质进编存中学纸文容保解留和论内法了保版许的方 全权予允，叉他 完有电并论其人校的制位或 本学文复学印 咩 峻日 畏2 章1 b _ - 4 - - 寸a n a v i e r - s t o k e s 方程是最有实际意义的偏微分方程之一，原因是它描述的是 一种非常普遍的现象流体的运动。它有非常广泛的应用，倒如： 汽车、飞机和宇宙飞船的工业设计； 血液在血管中流动的医学研完； 空气、河流和海洋中的环境污染分析； 海洋和气候运动的建模、计算机模拟和预测。 n a v i e r - s t o k e s 方程在具体的解中通常呈现极其大量的高频细节。更具 体地来说，流体的速度动态地、频繁地在很多地方出现陡变。极端的情况 用数学术语来说是阎断，比如激波( 想象一下我们在海滩上看到的一波叉 一波的海浪) 。n a v i e r - s t o k e s 方程最一般的形式是 p ( 基+ 剥警毒n 叭， 其中，对流加速度是这些陡变的主要原因，这是因为具有大r e y n o l d s 系数 的对流会在不同速度的流体的接触面处产生速度场的间断。更具体的论述 可以在引用1 6 0 j 中找到- n a v i e r - s t o k e s 方程的解申这些陡变的特征和规模是同时需要自适应网格 和并行计算的原因。这些陡变的特征是其动态随机性和尖锐性，这使得 1 静态网格很难适应其数值求解用均匀网格进行的数值求解，即使只是 中等精确程度的要求，都需要非常细密、未知量非常多的网格。非均匀 网格仅在这些陡变的分布已知的条件下有用，这种条件在数值解求出来之 前是很难达到的。因此，动态的自适应网格比静态网格更适合这种计算。 然而，光有自适应网格还不够。试想，这些陡变发生在那么多的地方而且 又那么频繁， 即使自适应网格已经大大减少了未知量的个数，这些大量 的高频细节还是会使网格很大的。因此，这些自适应网格计算还必需并行 化。 有限体积元方法( f v e m ) ，又称广义差分方法1 8 ，1 0 ；1 1 i1 5 ；1 6 , 2 4 ，2 6 ， 2 8 ，3 4 ，3 9 ，4 5 ，4 6 ，5 5 ，5 6 ，6 6 ，7 5 l ，是一种很适合n a v i e r - s t o k e s 方程数值求 解的计算方法。简单来说，有限体积元方法结合了有限差分法( f d m ) 和有限元法( f e m ) 的长处。它既有有限差分法处理守恒律和陡变的方便 性，又有和有限元法相类似的优美的分析理论和接近的几何灵活性。本论 文在稍后会简要地介绍一下这个计算方法。更多的介绍可以在引用1 4 5 l 中 找到。 因为上述原因，本论文的工作主要是有限体积元的自适应并行计算。 然而，对本论文而言，实现有限体积元的自适应并行计算仅仅是出发点。 本论文的中心，同时也是本论文的特殊之处，在于解决一个非常有时代 性问题，而且其解决方法是非常有创新性的具体地来说，这个极具时 代性的问题就是多核处理器出现后，现行最主要的并行计算方式集群 计算必须求变，否则无法适应多核处理器的发展和进步。目前文献 记载的两种方案( 见1 4 ) ，第一种方案需要全方位地改变硬件、操作系 统、集群计算作业调度系统和用户程序等多个层面，由于其工作量巨大， 这个方案在短期内很难实现第二种方案需要彻底地改变用户程序的算 法。由于有限体积元对网格质量有特殊要求，致使其算法很难按照第二种 方案的做法去重新设计，所以第二种方案对有限体积元的自适应并行计算 意义不大。既然这两种方案都不适合，本论文就提出自己的解决方案。 本论文既不改变硬件，也不改变用户程序，只重构并行计算系统运行库中 基本通信模块( 包括数据交换、并行求和以及并行求前缀和) ，并对用户 程序进行重新编译。经理论分析和实践验证，本论文的解决方案相当简 单实用，很适合于有限体积元的自适应并行计算。 多核集群计算是一种刚刚出现的计算方式，但它以后会成为集群计算 的主流。集群计算是在计算机集群上的并行计算。当前排行前5 0 0 的超级 计算机1 1 l 绝大部分都是计算机集群。所谓计算机集群，是由局域网连接 在一起的一组计算机，其中的每一台计算机都叫做这个集群时的一个节 点在过去，传统的计算机集群是每个节点配备一个单核的处理器。从 最近几年开始，随着多核处理器的出现，集群里的节点开始配备多核处 理器。这些多核处理器在一个芯片上嵌入了多个计算单元，其中的每个 2 计算单元都像过去的单核处理器一样功能齐全。 多核集群计算跟以前的传统的集群计算是有很大区别的。美国加州大 学伯克某分校在2 0 0 6 年( 1 4 1 ) 预言，传统的集群计算程序只可能在每节点 核数不超过8 个时得纠性能提升，超过每节点8 个核后就很难再提升了。 这一预言在2 0 0 8 年得到证实。在2 0 0 8 年，一篇i e e es p e c t r u mo n l i n e 的报 导( 1 5 0 1 ) 揭示， 美国国家实验的工程师们模拟了一些重要并行计算程序 在多核集群上的运行，发现其性能同时受制于集群网络和节点内存的性 能。也就是说，除非网络和内存的性能提高能跟得上处理器的发展，否 则这些程序很难在每节点核数提高的时候改进性能，甚至会反而变幔。 很不章，网络媒介不是半导体，在发展上不像处理器和内存那样受摩尔定 律管制，因此多核集群的网络问题会成为一个很严重的瓶颈。 1 1多核集群计算的困境 引用1 5 0 1 中所提到的那纽工程师模拟了未来的每节点8 、1 6 和3 2 核的计算机 集群，用来测试当前某些典型的并行计算应用程序在这些集群上的性能。 测试出来的结果相当令人沮丧。他们发现有限的网络和内存带宽使得这些 程序的性能并不随着每节点核数的提高而提高，甚至在数据访同量很大的 情况下反而降低。 简单来说，这个困境是由同一个节点上多个核共用有限的帝宽引起 的。如果一个节点有p 含核，这些核同时通过网络访问大量的数据，那 么每个核所得到的帝宽平均下来就不走过最大带宽的1 p 。更严。重的问题 是，假设集群里有n 个节点，每个节点上有p 个核，并且每一对不同的核 之间都有一个通信连接，那么每个节点就同时有p ( n p 一1 、个通信连接， 从而每个通信连接所能用纠的带宽平均下来就不超过1 ( p ( n p 一1 ) 1 最大带 咖 嚏 1 2处理器的多核发展是不可避免的趋势 既然多核集群计算有这样的困境，人们很自然会同：为什么处理器一定要 向多核发展呢? 有没有替代的方案，不用增加核数都可以提高处理器的性 能呢? 很不章，答案是否定的。原因是，处理器过去的发展方向已经碰 到了一些物理极限，使得现在不得不向多核发展。这些物理极限包括“内 存墙”、“指令级并行性墙”和“耗电量墙”。 过去，处理器发展的主流趋势是提高处理器的工作频率。例如，1 9 7 8 年 出现的8 0 8 6 处理器是5 到1 0 m h z ，1 9 8 2 年出现的8 0 1 8 6 处理器是6 刭1 2 m h z ， 3 1 9 8 9 年出现的8 0 4 8 6 处理器是1 6 刭1 0 0 m h z ，1 9 9 3 年出现的p e n t i u m 处理器 是6 0 刭3 0 0 m h z 2 0 0 0 年出现的p e n t i u m4 处理器1 3 到3 8 g h z ( 1 g h z = 1 0 0 0m i a z ) 直到2 l 世纪初，处理器发展的主流趋势才因为下面的技术 原因转变为在一个芯片上嵌入越来越多核。 “内存墙 在1 9 9 5 年就讨论到了，见w u l f m c k e e 的文章1 7 2 卜这个问题 是由处理器速度和主内存速度这同的差距造成的。由于处理器比主内存 要快得多，处理器的执行经常会被卡住，以等待从主内存过来的数据传 送。一个减少这种等待的方法是在处理器内部嵌入高速内存作为缓冲区， 把经常用的数据存在这个高速缓冲区里，从而不用每次都从主内存拿数 据。这个缓冲区叫作处理器的高速缓存。要注意，高速缓存通常比主内 存小很多，原因是造价过高。因此，处理器还是会有在高速缓存上找不 到需要的数据的时候，从而还得时不时到主内存拿一下数据。一个重要 的观察结果是，当处理嚣速度和主内存速度速度的差距继续拉大时，处 理器需要刭主内存拿数据的频率就会增加，除非高速缓存的大小也相应提 高。可是，高速缓存是如此的昂贵，这绝对是继续提高处理器工作频率 的一大阻碍。 “指令级并行性墙，t 是指w a l l 在1 9 9 1 年提刭的指令级的并行性的局$ 1 7 1 1 开采指令级的并行性是提高处理器性能的另一途径。典型的例子是现代处 理器中的超标量流水线处理。在那样的一个处理器里，每一条指令都被 分解为一个微操作序列，而这个序列是在处理器里的流水线式硬件里执 行的，从而不同的指令可以在执行时间上重合，达到某种程度上的并行 性。然而，指令级的并行性非常有限并且与应用程序相关，因此处理器 在这方面的性能提高这几年已经差不多挖刭尽头了。 “耗电量墙”与处理器日益提升的耗电量有关。2 0 0 1 年，i n t e l 的工程师 们在一篇文章( 1 3 3 1 ) 中提出，靠提高处理器的工作频率来提高性能是越 来越不节电了。2 0 0 5 年，h o f s t e e 的文章( 1 3 6 1 ) 揭示，靠指令级并行性来 提升性能同样是很浪费电能的。 所有的这些技术问题都强迫主要的处理器生产商改变处理器的发展策 略，转用芯片级多处理器( 也就是多核处理器) 的办法。这个办法可以 追溯到s t a n d f o r d 大学1 9 9 6 年的一篇文章( 1 5 1 1 ) 人们认为这个办法可以 有效地减轻这述技术问题所带来的影响，尽管它要求计算机多层次盼变 革，特别是在软件层次上。 1 3多核集群计算困境的解决刻不容缓 既然多核发展会统治未来的处理器发展趋势，我们应该在什么时候想办法 解决多核计算的困境呢? 多核发展的预测暗示，这一困境的解决是越快 4 越好，因为处理嚣里的核的个数会增加得非常快。国际半导体发展规划 ( 1 t r s1 3 8 i ) 2 0 0 7 年的报告预测，处理器的核数将以每两到三年翻一翻曲 速度增长。s t e n s t r o m 在他2 0 0 6 的文章1 6 2 i 里说，处理器的核数是每三年翻 一翻。美国加州大学伯免莱分校则预言( 1 4 1 ) ，处理嚣的核数是每个半 导体换代翻一翻。基于这些预测，p a r a 0 6 会议的一个w o r k s h o p 定论，为 多核计算机重新研制计算数学软件的时候已经刭了。类似地，a n w a r ， i n t e l 微处理器的主工程师，在2 0 0 8 年六月3 0 日建议1 3 1 1 人们马上为几十个 核、几百个核甚至几千个核的处理器设计并行程序，要赶在这些处理器来 到之前做好准备。 1 4现有文献中可查刭的相关解决方案 现有文献中可以找列两类比较完整的解决方案。还有一些只解决其中一些 很具体的相对较小的问题的技术，这里不一一列出了。 w a r f t ( w a r e nr e s e a r c hf o u n d a t i o n ) 基金推动一项在印度进行的 研究，以处理多核集群计算中的困境。这项研究在2 0 0 7 年发表了两篇文 章1 6 9 ，7 0 l 。第一篇文章提议将主内存移入处理器，从而节点上的内存带 宽多核争用紧张度可得到缓解。第二篇文章建议用一个全新的执行模式来 执行集群计算，其中需要多层次的重构，包括操作系统、计算作业调度系 统和计算程序的全方位重新设计，同时在所有节点上运行多个计算作业， 从而在网络层次缓解多核争用紧张度。这两篇文章的方案复杂度相当高， 工作量相当大，可行性有待考证。 v a l i a n t ，b s p 并行计算模型的发明者1 6 7 1 ，将b s p 模型修改以适应新的 多核集群计算。并行计算模型是用来设计并行算法的， 主要作用是提供 统一的标准来衡量并行算法的复杂度，方便算法的设计。他的新并行计 算模型称为m u l t i b s p ，具体内容可以在他的一篇2 0 0 8 的文章1 6 8 l 找到。这 个新模型的基本思想是：新出现的多核集群的内存结构是多层次的， 包括 网络( 全局) 层次、节点层次、核层次( 指高速缓存) ，而且每个层次可 以在固定时间里存取的内存的大小都是有限制的， 因而算法在设计时亦应 该相应地多层次并且符合每个层次的参数和条件。这个方案也许对很多问 题是适用的，但： 重新设计算法是不容易的，更何况多层次的复杂度分析相当复杂。 给定一个问题，例如下面即将介绍的有限体积元自适应计算，是否一 定存在符合每个层次的参数和条件的多层次算法呢? 这是个疑问。 5 1 5 有限体积元的简介 有限体积元方法可以看作是一种使用双重网格的广义p e t r o v - g a l e r k i n 方 法。 单元网格用于试探函数空同，用分片多项式来逼近数值解。 控常j 体积网格用于检验函数空闻，这个空间包括了控制体积的特征函 数。 1 5 1网格的双重性假设 简单而言，如果下面的条件成立，一个单元网格和一个控制体积网格就构 成一对双重网格：给定标准单元和将标准单元映射为网格单元的仿射变 换，每个网格单元中的控制体积分划( 也就是控制体积网格与这个单元的 交) 相同于标准单元上预定义的控制体积分戈i j 在相应仿射变换下的映像。 下面是更详细的解释。 假设区域q 被划分为单元网格丁= f k ，其中k 是网格单元。证= 毪) 表示未知变量的节点集合证丁+ = 蟛 表示控制体积网格， 其 中雕是包含猕的控制体积。让詹表示标准- t - 元，并且证对，争表示 在露的节点集合和控制体积网格。如果下面的条件成立，则丁和r 是双重 网格：所有这些仿射变换 r、 厅：k t ，厅( k ) = k lj 都满足关系 厅( m = kn ， 厅( r ) = k np 注意：每个网格单元都有多个厅的选择，这取决于露的哪个顶点被 映射到k 的哪个顶点。上述假设应被理解为与这些选择无关。 1 5 2用有限体积元方法离散化方程的一个例子 为了讨论的方便，请考虑如下的模型问题。下面这个问题女 n a v i e r - s t o k e s 方程中的扩散( 粘性) 项有密切关系。因为扩散项的长距离效 应，n a v i e r - s t o k e s 方程中只有扩散项需要本论文讨论的与网格质量相关的 特殊处理， - a u ( x ) = f ( x ) v x q ， v u ( x ) n ( x ) = 0 y x 凇 6 ( 1 5 1 ) ( 1 5 2 ) 违里，如囤l5 1 所示，n 是r 2 中以扰2 为边界的多过形区域 f 是一个光滑 函数，n 是指向外面的法向量场。遣一方程的解在下面这个商空问里存在 并且唯一： w 【等价于w 2 昔wj w 2 是常量 设标准单元为单位三角彤，开即 霄= ( z ，p ) r 2 ：。兰0 ，0 $ + p 茎1 让对表示膏邵三个顶t 曲集台，并且设疗被划分为控制网格于，如图l5l 的 古上角所示。则f l 上的双重网格的生成如同一个圜中的左下角所示。 f i g u r e15 1 。左上：单元网格- 右上：被划分的标准单元。左下用标准 单元生成7 控制体柱网格的双重网格。 在运对盈重网椿上- 取试探函数空间为由标准单元上的线性函数在边 界条件下生成的基函数空问，而检验函数空间m l 取为控制体积上的常量函 数。这样， 方程( 1 5 1 ) 和( 1 5 2 ) 就可以被离散化为寻求试探函数空间中 的u 丁，使得 n 丁( 札丁，v t ) = ( ，秒丁) 对所有检验函数空间中的t ，p 成立，其中 并且 。丁( 岍，协) = a k ( w t ，即) ( 1 5 3 ) k 丁 a ( 【w t ，v t j = 埘p 0 v w ：咖非常量 j 本论文称这一条件为网格的全局代数约柬。 8 这个全局代数约束是很难在实际中应用的，除非它被分解为一组单元 上的约束。事实上，这种分解是可以做到的。全局代数约束的一个充分 条件是用下面这个由单元分析所导出的条件：对任意网格单元k 丁： w ( a 乏+ a k ) w 0 v w ： q 奶在k 上非常量， ( 1 6 1 ) j 其中a 耳是( 1 5 4 ) 中的o k 的离散化形式： a k = f a k ( 咖，讥) 】 这些只需要察看单个单元的a 在分析上明显易于a 丁。 在程序实现上，后面即将介绍刭的正曼l j 自适应算法对给定的有限体积 元方法有一个可接受网格质量的下界，而这个下界用单元约束来描述比用 全局代数约束来描述更为合理。初看时可能会觉得条件( 1 ：6 1 ) 是全局代数 约束的充分非必要条件。然而如果动态地来看网格自适应的过程就不同 了，这是因为在给定初始网格之后，这个自适应算法只使用有限个单元 形状来动态地局部加细或减细网格，而初始网格中的每个单元的形状都有 可能因为局部细分而加大它对全局的影响，甚至会起决定性作用。因此 下面这个定理说出了一个事实：如果要通过对初始网格的控制来保证动态 网格的质量，那么( 1 6 1 ) 中的单元约束就是必要条件。 定理1 6 1 风格单元k 的形状泱定a k ，但其大小与a 无关。 这个定理实际上是后面将要介绍的定理2 1 6 的直接结论。 1 6 2 有关的工作及结论 网格质量问题其实是与有限体积的高次元格式紧密联系的。高次元格式使 有限体积元方法可以用较少的未知量个数达到较高的精度，这对高精度 的n a v i e r - s t o k e s 方程的求解来说是很重要的。同时，后面我们会看到，高 次元的高阶格式是需要本论文所讨论的这种特殊的网格质量控制的原因。 早期有限体积法高次元的研究包括c h e n 的文章1 1 9 ，1 8 l ，p h e m o n o t 的论 文1 2 5 1 ，l i e b a u 的研究1 4 7 1 ，l i ，c h e n 和w u 的著作1 4 5 1 ，p l e x o u s a k i s 和z o u r a r i s 的 工作1 5 3 l ，c a i ，d o u g l a s 和p a r k 在矩形单元上的研究1 1 4 1 ，还有x u 和z o u 的文 章1 7 4 1 本论文在网格质量方面的工作受弓l 用1 4 5 1 所启发。该引用还给出了二 次、三次元的网格约束的一些充分条件。另外，引用1 7 4 1 也给出了若干二 次元的约束条件。本论文与这些约束条件的最大不同之处在于，本论文 力求发挖掘这些约束条件底下的运算规律，从而不需要人为的分析，都 9 可以把这个约束条件算出来。这是基于程序实现的角度出发的。因为有 一个算法可以处理一大类的计算格式，就不需要单独对每个不同的计算格 式分析和编程了 1 7自适应计算的简介 典型的自适应计算迭代地执行下面的步骤： 1 在当前的网格上计算求解。 2 做一个后验估计。如果误差已经足够小，则算法结束；否则根据估计 所得的误差选取并且标记要加细或者减细的网格单元。 3 根据上一步骤的标记动态地改变网格，在某些地方加细，在某些地方 减细 4 跳转到第一个步骤。 为了保证在自适应的过程中的网格质量，一个特殊的网格自适应技术 被采用；这一技术就是在单纯形网格( 特别是在二维的三角网格) 中广泛 使用的正则自适应法( 叉称红绿自适应法) 。重要的是，这一自适应法所 生成的网格单元都与原有的初始网格单元相类似，因而网格质量可以得到 保证。 1 7 1二维情况下的正则自适应法 正则自适应算法是由b a n k ，s h e r m a n 和w e i s e r 在1 9 8 3 年提出的1 3 l ，并在一 个名为p l t m g 的软件包中成功地得到应用1 5 l 。在随后的一些文章中这个 算法得到推广1 6 ，9 ，1 7 ，4 8 ，7 6 l o 在二维在情况里，如图1 7 1 所示，被标 记要加细的三角形被连接三边的中点，从而正更i j 地分为四个子三角形。 为了使得刭的网格保持协调，如图1 7 2 所示， 与这些三角形( 图中红色 者) 邻近的三角形( 图中绿色者) 要按照相邻的边来二分为两个子三角 形。在原来的算法中，这些二分三角形的子三角形不能继续细分，除非 首先得把二分的三角形转为四分。这个原来的算法会遇到一个比较复杂的 情况一个没被正昊f j 细分的三角形是有可能同时与两个被正则细分的三 角形邻近的， 因而这个三角形在细分的过程中必定会先被二分，再被四 分这在串行计算中问题不大，但在并行计算中很难高效实现。幸好， c a r t e n s e n 在弓l 用1 1 7 l 中提出对这个算法的修改。在同时相邻两个正昊i j 细分 三角彤的时候，二分三角形的其中一个子三角形可以被继续二分一次。 1 0 在后面的算法中奉论文特使用c a r l e l s o n 修改过的这个算法。并且使用 并行的柑数据结构来* 现它。后面鲁有详细介绍。 f 墉r e17 i ：三角形的正劓细分 影勺y 、 f i g u r e l72 ：与正州细分三角彤相邻的三角形要二分来保持网格的协调性 1 7 2 其它相关的罔格自逢应法 1 9 8 4 年，r i 提出三角形网格的最长遗二分自适应算法5 4 i 。这个算 法细分三角彤时，走；的是选择三角夥的最长边来二分。有小数三 角形要选短一点的遗来二分，“确保网梏的协调性。这个算洼也能保 1 正网格的质量丰论互之所以恐有选用这个算法是因为它的并行宴 现i j ( ) i 比c a r l m h e n 修改过的正列自适应珐的并行窭现多t 一些中间步骤 ( 比如选取非依鞋集) ，性能上要慢一些。另外，后面我们套看到有 限体积元的同格质量控制用三角彤的边长之比垂是活，而最长迪二分法的 最差情咒估计是甩最一1 、角来做的1 5 2 1 相对而言迁那么适合有限俸积元计 算。 1 9 7 9 年，s e w e l 提出三角彤r q 格中的最新节点二奇自适应算法5 。1 。在 返十算法里，三角形的二争是拄照它的最新节点来做的。这个算击在实 际丘用中有些清况很难处理。例如算法中要求细分的三角彤必须是一对 一对的，从而有可能有些无法配对的三角彤得不到细分m i t c h e l l 4 9 f 中 去掉7 这个要求，但不幸的是这十修改导致那些无法配对的三角彤必颈递 归地进行细分。运在并行计算中很难娈现。 1 8基本通信模块重构的简介 重构数据交换、并行求和、并行求前缀这三个基本通信模块和是本论文解 泱多核集群计算的方法。选取这个方法的原因主要有两个方面。首先， 这些基本通信模块是并行计算中通信的重要形式。其次，重构这些基本 通信模块能够以非常小的代价解决多核集群计算的网络负荷困难。 1 8 1基本通信模块的重要地位 数据交换、并行求和、并行求前缀和是并行计算中通信的三种主要形式。 其它形式包括同步f 1 3 l 、任务转移1 5 7 ，6 4 l 和主从调度7 l 。这三种基本通信 模块和同步在使用经典并行编程模型的应用程序中很常见。经典的并行编 程模型包括p r a m 2 7 ，4 2 ，6 3 ，b s p 3 5 ，6 7 1 和l o g p 2 2 1 。然而， 因为在很多 并行编程a p i 中的标准化和优化( 在m p i l 3 2 ，6 1 l 和p v m 2 9 ，3 0 1 中) ， 同步 比数据交换等基本通信模块受多核处理器的影响要小得多。其它的通信方 式比较少见，只在像分子动力学模拟那样的应用程序中常用。 1 8 2数据交换重构的基本思想 本论文重构数据交换的基本思想有三方面： 尽可能地减少通信连接，从而抑制多核处理器对集群计算的负面影 响； 把通信连接分配给多些核一起来做，从而增强多核处理器对集群计算 的正面作用； 用不同的核处理纠不同节点的通信连接，从而第一方面和第二方面可 以无冲突地实现。 1 8 3并行求和、求前缀和重构的基本思想 这两个基本通信模块的重构可以分三步进行： 优化本地计算的协调过程； 优化全局通信的性能： 处理节点上网络共用紧张度的问题，并且更大程度地利用本地的共享 内存。 1 2 1 8 4相关的工作 用减少通信连接来抑制多核处理器对集群计算的受面影响的主意来自于弓l 用1 4 1 l 。把通信连接分配给多些核一起来做来增强多核处理嚣对集群计算 的正面作用的主意来源于引用1 4 3 l 。然而，在实现方法上，上述两个引用 提出的方法是冲突的。第一个引用的方法在每个节点上只用一个核负责对 外通信。第二个引用的方法用节点上的所有核共同对外通信。本论文的 方法在重构数据交换运一部分中的最大特色是用_ 个全新的任务调度，来 解决这两方面主意表面上的冲突。 1 9正则自适应的树结构并行算法简介 1 9 1正曼l j 自适应网格的树表示 在本论文的并行实现中，自适应网格是存放在若干层共享数组里面的。 最上面一层数组存放的是初始网格的三角单元。第二层存放的是初始单元 经过一次正则细分后产生的子单元。每三层存放的是第二层的单元进一步 正则细分后产生的子单元。如比类推。 每个三角单元在数组中都占用一个元素的内存。这个元素里面包含了 相应单元的父单元和邻近单元的序号( 位置) 。除了把指针换成序号， 整个数据结构就像串行计算中树结构的一样。 1 9 2自适应过程中的基本操作 本论文的树算法的基本操作实际上只有两个：( 1 ) 并行细分算法：对已 经标记了要细分的单元进行细分，并且保持其网格协调性。( 2 ) 从网格 生成方程：并行地遍历这个树结构，把这个树结构平面化戍一个数组，从 而进一步产生离散亿的方程。 1 9 3相关的工作 本文的树遍历算法是基i f - r i c h a r dc o l e 和u z iv i s h k i n 的并行树处理1 2 1 l 而设 计的。基本思想是，科用欧拉游路( e u l e rt o u r ) 技术把树的遍历转化为 链表算序问题，而链表算序问题又可以避一步用指针双跳算法变为简单的 几个数据交换加上加法和赋