时域波动方程边界元法的优化与并行化

时域波动方程边界元法的优化与并行化引用此版本：Berenger Bramas.时域波动方程边界元法的优化与并行化。分布式、并行和集群计算。波尔多大学，2016年。英语NNT：2016BORD 0022。电话：01306571HAL Id：tel-01306571https://theses.hal.science/tel-013065712016年4月25日提交HAL是一个多学科的开放获取档案馆，用于存放和传播科学研究文件，无论它们是否已这些文件可能来自法国或国外的教学和研究机构，或来自公共或私人研究中心。L’archive ouverte pluridisciplinaireTHAPOSE法国总统府L’UNIVERSITÉ DE BORDEAUX数学与ParBRAMAS BérengerPOUR OBTENIR LE GRADE DE DEDOCTEUR专业：信息Temporel地区边界要素平衡方法的优化与改进搜索：15 Février 2016报告员抵达后：BIROS乔治·埃塞尔-得克萨斯大学奥斯汀分校VUDUC理查德·埃塞尔（Associé）-佐治亚理工学院在委员会对以下作曲家进行审查时BIROS George Cherseur-德克萨斯大学奥斯汀分校特别报告员COULAUD Olivier Directeur deRecherche - Inria Bordeaux-Sud-Ouest Directeur de ThèseTERRASSE Isabelle Directrice de Recherche -空中客车集团审查员VUDUC Richard Rupseur（Associé）-佐治亚理工学院特别报告员2016时域波动方程边界元法的优化与离散化贝伦热·布拉马斯在履行的要求博士学位的主题计算机科学二月15th，2016裁判乔治·比罗斯·理查德·武杜克教授-德克萨斯大学奥斯汀分校副教授-佐治亚理工学院委员会乔治·比罗斯·奥利维耶·库洛·帕斯卡尔·哈韦·斯特凡·兰泰里·雷蒙德·纳梅斯特·纪尧姆·西尔万·伊莎贝尔·特拉斯·理查德·武杜克教授-德克萨斯大学奥斯汀分校研究主任-InriaBordeaux-Sud-Ouest研究员- IFP Energiesnouvelles研究总监-Inria Sophia Antipolis教授-波尔多大学研究员-空中客车集团创新研究总监-空中客车集团佐治亚理工学院（Georgia Instituteof Technology）审查员委员会主席审查员报告2©2016版权所有数学与计算机科学博士研究生院计算机科学系，Inria Bordeaux-Sud-Ouest，Université de Bordeaux1.已于2016年02月18日星期四10：47IIITemporel地区边界测量方法的优化与改进贝伦热·布拉马斯简历边界元法（BEM）用于声学和电磁学中，以模拟具有时间分辨（TD）的声波传播她同意为多位朋友提供一个解决方案在此基础上，我们研究了不同角度下模拟TD-BEM的Nous décrivons le contexte denotre étude et la formulation utilisée qui该系统是利用运算矩阵/矢量creux（SpMV）进行自然计算的我们的工作围绕着SpMV的局限性，研究矩阵的排列，以及CPU上的向量化和GPU上的块的实现。 Nousmontrons que cet opérateur n'est pas approprié pour notre problème et nous proposons dechanger l'ordre de calcul afin d'obtenir une matrice avec une structure particulière.这一新结构被称为一个矩阵部分，并由一个专门的操作员进行计算我们对高性能计算的现代建筑的实施进行了该模拟结果与一种混合的方法（备忘录/分发）并行，以新的方法为基础，以均衡处理者之间的工作。这是一个复杂的二次曲线，我们可以通过多极快速方法（FMM）来研究它的加速度。我们利用不同的范式对FMM算法的并行化进行了大量的研究，并对FMM算法的执行机制进行了评述，以揭示FMM的潜力 恩，我们提出了由FMM加速的模拟TD-BEM的初步结果。Mots-clés ： calculhauteperformance ， programmationparallèle ， optimization ，vectorisation ， GPU ， methodedesélémentsfrontières ， équation desondes ，acoustique，éelectromagnetisme.III时域波动方程边界元法的优化与离散化贝伦热·布拉马斯摘要本文用声学和电磁学中波动方程的时域边界元法模拟波在时间上的传播它允许一次求解得到多个频域结果在这篇论文中，我们研究了一个有效的TD-BEM求解器的实现使用不同的方法。我们描述了我们的研究背景和TD-BEM制定表示为一个稀疏的线性系统组成的多个相互作用/卷积矩阵。这个系统自然是使用稀疏矩阵向量积（SpMV）计算的 我们通过查看矩阵重新排序以及在CPU上使用矢量化（SIMD）和在Nvidia GPU上使用高级块布局的SpMV内核的行为来研究SpMV内核的限制。 我们表明，该算子是不适合我们的问题，然后我们建议重新排序的原始计算得到一个特殊的矩阵结构。这种新结构称为切片矩阵，使用自定义矩阵/向量乘积运算符进行计算我们提出了一个优化的CPU和Nvidia GPU，我们描述了先进的阻塞方案上实现该操作符由此产生的求解器是并行的异构节点以上的混合策略，并依赖于一个新的启发式处理单元之间的工作平衡。 由于这种矩阵方法的二次复杂性，我们研究了使用快速多极子方法（FMM）的时域边界元求解器。 我们调查的并行化的一般FMM算法使用几个范例在共享和分布式内存，我们解释了现代的运行时系统是如何非常适合表达的FMM计算。最后，我们调查的实施和具体到我们的TD-BEM的FMM内核的参数化，我们提供了初步的结果。关键词：高性能计算，并行编程，优化，矢量化，GPU，边界元法，波动方程，声学，电磁学。IIITemporel地区边界测量方法的优化与改进贝伦热·布拉马斯艾滕杜简历利用声电磁领域边界条件的方法模拟时域边界条件的传播。她希望通过一个简单的解决方案来获得多个成员的结果在此基础上，我们研究了不同角度下模拟TD-BEM的Nous décrivons le contexte de notreétude et la formulation utilisée qui peutan=（ M0）−1（ln−Kmaxk=1Mk·an−k）.（一）该系统是自然计算的，它利用的是产生的矢量矩阵（SpMV），KMaxmatricesMk.在这方面的研究，我们intéressons在一个首要的时间aux limites du SpMV，de la置换矩阵在一个工业模拟器中，置换算法的计算结果更能与新的SpMV的增益进行比较我们的montrons qu 此外，这些块是由超过50%的零组成的，而SpMV的性能也是如此。 这一观察是解决我们问题的主要动机。在最初的公式中，外部边界是在KMax矩阵上，内部边界是在Chaque矩阵的线和列上。我们提出的réordonner的计算是为了获得矩阵的结构，特别是由一个非空的向量组成的矩阵;外部的边界在矩阵的列上产生在获得一份适当的备忘录后，这些矩阵不是由一种新的称为矩阵的结构所取代的Nous étudions总之，这些操作员对优化非常明智Nous décrivons notre implémentation sur CPU utilisant la vectorisation enC mais aussi en browsavec pour objectif de réduire le nombre∑iv记录回忆录。 我们设计了一种混合并行算法，用于我们的模拟，并证明了M0偏差批判的分解/求解是在增加nøud的名称。Nousétendonsnotresimulateurpourfonctionnersurlesarchitectureshybridescomposées de plusieurs CPU/GPU.Le calcul des matrices en transche sur GPU estexigeant，et il est obligatoire 我们比较了两种块结构和GPU上的Noyaux结构，并介绍了一种新的算法，以便在初始迭代之前使各单元的计算工作更加均衡。 Nous validons notre simulateur sur 5noeuds chacun composés de 24 CPU et 4 GPU，et nous démontrons que nous avons unebonne accélération lorsque le problème passe en memoire dans les GPU.Néanmoins ，cette approche matricielle a une complexité quadratique pour la résolutionmais aussi pour la construction des matrices. 这种复杂性可能存在局限性，因为我们可以利用多极加速方法来研究模拟的TD-BEM。Nous nous intéressons dans un premier temps a La parallélisation de la FMM est simpleavec une approche fork-join comme celle proposée par OpenMP. 在另一个角落，我们的蒙特利尔人评论说，执行机制非常适合于更好地探索FMM的发展。此外，我们的montrons评论améliquiles performancesgrâce à la propriété commutative des opérations de la FMM。我们建议采用一种新的结构，即采用一种新的执行机制，使FMM并行化 Cet arbre-en-groupe permet de parametrer la granularitédes tâches et potentiellement de réduire le nombre de dépen-dances. 我们介绍了一种新型的FMM，它以我们的组轴为基础，以记忆为基础。Nous comparons les comporcidence de différentesimplémentationsenmemoirepartagéeetdistribuéepourdessimulationsdeparticules.ScalFMM是FMM的一个自由软件包，在ScalFMM的书目中包含了不同的逻辑描述。然后，我们学习如何使用FMM和适当的FMM来实现TD-BEM求解器我们利用矩阵方法和相互作用矩阵来计算各层之间的相互作用，并利用FMM来计算各层之间的相互作用我们可以不使用所有的交互作用矩阵然而，FMM可能会采用不同的方法，并在参数上加以纠正FMM的不同运算可以在时域或频域计算，这需要傅立叶变换的补充。Nous montrons qu 此外，我们不需要计算整个FMM，因为它可以在计算时间时进行迭代。v内容介绍11 问题陈述31.1工业应用环境。 . . . . . . . . . . . . . . . . . . . . . . . . . . . .31.2波动方程的时域边界元法（TD-BEM）。. . .51.2.1线性系统公式. . . . . . . . . . . . . . . . . . . . . . . . .51.2.2相互作用/卷积矩阵。. . . . . . . . . . . . . . . . . . . . . .61.2.3分辨率算法。. . . . . . . . . . . . . . . . . . . . . . . . . . ..71.2.4研究测试用例。. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .81.3现代HPC。. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .91.3.1并行和高性能计算。. . . . . . . . . . . . . . . . .91.3.2科学计算中的数值运算112 背景和技术现状132.1 稀疏矩阵向量积（SpMV）132.1.1稀疏矩阵存储调查142.1.2重新排序202.1.3并行SpMV212.1.4加速器上的SpMV222.1.5实验实施例222.1.6摘要.232.2 并行计算系统242.2.1来自数据流的DAG的表达式242.2.2其他数据访问模式242.2.3StarPU252.3 快速多极子方法（FMM）272.3.1N体问题和直接计算272.3.2分层方法282.3.3快速多极子方法（FMM）282.3.4空间填充曲线302.3.5ScalFMM342.3.6其他并行FMM35vi2.4 波动方程的其他TD-BEM公式362.5 捐款373 TD-BEM矩阵方法393.1 SpMV-39型的优化尝试3.1.1矩阵排序和未知数排序393.1.2重新排序403.1.3重新排序一组矩阵443.1.4SpMV实现453.1.5未对齐块坐标存储（UBCOO）453.1.6UBCOO SpMV45的CPU实现3.1.7CBZ SpMV47的GPU实现3.1.8SpMV使用总结483.2 重新排序求和计算493.2.1订购邮政493.2.2切片属性503.2.3多步骤切片计算算法513.3 CPU52上的多向量/向量积3.3.1SIMD多矢量/矢量积533.3.2内存和组件优化543.3.3管理行向量的变化553.4 GPU56上的实现3.4.1GPU56上的切片计算3.4.2全面封锁做法573.4.3连续阻挡法593.5 企业化603.5.1同质节点的并行化策略603.5.2异构节点的并行化策略623.5.3并行线性解算器注意事项643.5.4总和的除法（远场近场）653.6 构建互动矩阵653.7 性能和数值研究663.7.1实验装置663.7.2平衡质量研究673.7.3时序触发率683.7.4测试用例723.7.5M073的线性解算器vii3.7.6平行研究753.7.7核心外执行793.8 矩阵计算总结804 并行快速多极子方法4.1 序贯快速多极子方法824.1.1算法824.2 共享存储器并行化834.2.1并行834.2.2tasks-and-wait任务和等待854.2.3第86节任务和等待4.3 混合并行化（MPI/OpenMP）874.3.1分布式内存（Full-MPI）874.3.2通信隐藏策略914.4 Tasks-and-mandencies任务和职责的简化924.4.1FMM直接无环图（DAG）934.4.2FMM94任务和职责4.4.3组树数据结构964.4.4基于任务的群组树994.5 分布式任务和权限1004.6 在基于运行时的FMM101中启用加速器4.7粒子模拟并行研究1014.7.1共享存储器并行化1014.7.2分布式内存并行化1044.8 摘要1045 使用FMM109的时域边界元求解器5.1 时域BEMFMM公式1095.1.1直接方法的局限性1095.1.2从过去接受还是从未来接受1105.1.3通过FMM八叉树对网格的空间划分1115.1.4原则和公式1125.1.5单位球离散化1165.1.6APS（近似长球形）1175.2 运营商1185.2.1P2MM2M120M2L123viii5.2.4L2L1245.2.5L2P5.3 优化1255.3.1关于时间或频率算子选择的讨论5.3.2不完整的FMM1265.3.3ScalFMM作为并行化引擎1265.44DFMM127的表达5.5初步数值结果和平行研究1275.5.1单位球面的表示法1275.5.2时移评估1285.5.3参数化1295.5.4测试用例1305.5.5时域与频域运营商........................................................................................................1315.5.6矩阵方法与FMM方法................................................................................................1325.6总结与展望137结论139前景141附录143附录A声学TD-BEM公式144附录B数学151附录C硬件概述160附录D FMM172附录E HPC软件工程182图189表格一览表190算法列表191源代码192参考文献202ix我把这篇论文献给我的妻子多萝西我很快奥兰多他们的支持和他们的耐心。我也把这篇论文献给我的家人，特别是我的父母，他们把我推向了这条道路。1介绍波动方程的时域边界元法（TD-BEM）用于模拟声波或电磁波在表面网格上的传播更准确地说，它模拟了源发出的入射波如何随着时间的推移在网格的离散化元素上传播和反射时域结果与傅里叶变换的转换允许我们研究几个频率，这使得TD-BEM求解器适合于研究宽带应用。然而，频域边界元法（FD-BEM）比时域公式[1]得到了更多的研究离散化元素的大小与研究的频率成比例地减小，因此覆盖网格的未知数的数量成二次方地增加从内存和工作负载的角度来看，使用超级计算机和高性能计算（HPC）是计算大型模拟所必需的。本研究中使用的TD-BEM公式已在[2]中引入，其中问题表示为具有使用稀疏矩阵向量乘积（SpMV）自然计算的求和阶段的线性系统。然而，这种方法有两个主要局限性。首先，众所周知，SpMV是一种内存受限的操作，它在大多数处理单元上实现了一小部分峰值其次，矩阵的生成和系统的求解相对于网格的离散化元素的数量具有二次复杂性，这可能限制我们能够计算的问题的大小在这篇论文中，我们研究如何解决这两个问题，并提供一个有效的TD-BEM求解器的现代HPC架构。由此产生的应用程序旨在取代现有的TD-BEM软件的基础上的SpMV。研究范围处 理 单 元 的硬件结构使得高效SpMV的开发极具挑战性，因为操作数量与数据大小的比率很低。在我们的例子中，通过重新排序计算，我们得到了每行一个稠密向量的特殊矩阵，这增加了这个比率。因此，我们专注于为CPU和GPU上的此类矩阵开发一个高效的内核由此产生的操作符不是一般的SpMV，而是特定于我们的问题。然而，硬件特性的开发和管理是经典的，可以应用于各种其他开发。在[3]中提出了一种快速多极子方法（FMM）内核来加速我们的TD-BEM。 我们研究了一般FMM算法的不同并行化方案，并将它们用于TD-BEM FMM内核的实现。 这个内核已经实现了高级别的优化，但它仍然是我们的TD-BEM的一个特定公式。曼尼托邦组织。第一章详细介绍了我们的声学TD-BEM公式，并描述了当今的HPC架构及其相关的挑战。第二章是相关工作的展开，它用它的2优化实现和一般FMM算法。最后三章描述了我们的贡献，目前的数值和性能研究。在第三章中，我们重点讨论了TD-BEM的矩阵方法然后，在第四章中，我们研究了FMM在共享内存和分布式内存中的并行化在最后一章中，我们重点研究了有限矩法加速的TD-BEM，并给出了初步的结果。最后，在附件中，我们提供了与我们不同实现的底层数学和算法相关的材料鸣谢：这项工作得到了空中客车集团、Inria和Conseil Régional d'Aquitaine倡议的支持31问题陈述在本章中，我们将概述与我们的工作相关的不同领域的背景和背景我们首先回顾一下，我们的开发是在工业背景下完成的，我们描述了波传播模拟的目的然后，我们给出了我们的时域边界元法（TD-BEM）的核心公式表示与稀疏矩阵矢量积（SpMV）。最后，我们总结了当今HPC的问题，解释了当前HPC应用程序开发中面临的挑战。1.1工业应用环境。这项工作是与空中客车集团创新公司合作完成的，空中客车集团创新公司是空中客车集团的一个实体，致力于为空中客车集团各部门（空中客车民用飞机，空中客车防务航天，空中客车直升机）的使用进行研究和开发。20多年来，数值分析团队一直致力于波浪传播模拟的积分方程和边界元方法。由此产生的软件解决方案是日常基础上的声学安装效果计算，航空声学模拟（在耦合方案与其他工具），并在电磁天线选址，电磁兼容性或隐形。自2000年以来，这些频域边界元法（BEM）工具已经扩展了多极算法（称为快速多极方法），使我们能够在合理的时间内在并行机上解决非常大的问题，具有数千最近，H矩阵技术已经使得快速直接求解器的设计成为可能，该快速直接求解器能够以非常高的精度解决具有数百万未知量的问题，而没有与迭代求解器相关联的通常缺点：不能控制迭代次数，难以找到好的预处理器等。与此同时，空中客车集团创新公司致力于时域边界元法（TD-BEM）的设计和优化本研究中提出的应用程序是工业计算工作流的一个层，必须足够强大才能在工业上使用我们将网格的管理或数据的生成委托给一些现有的黑盒来集中处理4他们的计算。 在算法或外部库方面的不同选择受到工业最终用途的限制。波动方程应用的边界元法波在声学和电磁学中传播的研究有许多应用。在航空工程中，它们用于设计飞机，但也用于保证其鲁棒性。例如，发动机的位置很重要，因为它会改变驾驶舱内听到的声音的音量此外，飞机的驾驶舱主要由金属组成，所有的通信都依赖于电磁信号，因此确保发送到发动机/机翼控制器的消息或与地面的通信是正确的至关重要必须避免组件之间的干扰，并且必须考虑到对闪电等外部源的反应。正如我们在下一节中所阐明的，为了研究小面/元素上的问题，对象被离散化。在图1.1中，我们展示了一个有三个天线的飞机网格示例图1.1：在网格中离散化飞机在图1.2a和图1.2b中，我们展示了电磁模拟的两个结果。第一个是用来研究天线对机舱的影响，第二个显示了一个模拟和实际措施所获得的结果之间的比较。(a) 从模拟结果来研究天线的发射(b) 通过模拟或与实际测量数值获得的波发射的比较图1.2：电磁学研究示例。5∑k=1∑1.2波动方程的时域边界元法许多物理过程可以用偏微分方程（PDE）来描述，但只有少数特殊情况可以解析求解为了克服这一缺点，许多数值方法已经发展起来，如边界元法，有限元法或谱方法。边界元法（BEM）是一种数值计算方法，当这些偏微分方程被公式化为积分方程（即，以边界积分方程形式-BIE）。这些方法是针对空间域上的边值问题转化为物理域边界上的积分方程所产生的边界积分方程而发展的在边界元法中，问题的三维几何形状在空间中被二维表面网格单元离散化，以仅描述原始对象的边界[4]。我们的问题公式的数学细节在附录A中给出，但它不是理解我们的贡献和实现部分的先决条件对物理方面和基本数学公式感到关注的读者可以阅读导致线性系统表达式的完整描述给出的问题描述和公式取自[3]，而时域公式最初是在[2]中为电磁学引入的电磁学中麦克斯韦方程组的线性系统与这里提出的声学系统相似因 此 ，在这项研究中提出的计算方法可以很容易地适应电磁学。1.2.1线性系统公式化在这一节中，我们回顾了我们的TD-BEM公式的关键点，并重点介绍了其线性系统的表达。 具有速度c和波长λ的入射声平面波w在边界Ω上发射。曲面Ω通过经典有限元方法离散化，该方法使用三角形单元，其中未知数/自由度位于顶点处我们用N表示系统中未知数的总数波动方程也在时间上离散化，步长为Δt，时间步长T的数量必须从研究的频率范围和所考虑对象的大小中选择 照明矢量ln包含入射波w在迭代时间t n= nΔt对网格未知数的影响。 一旦波照亮了定义未知数的位置，它就会被它们反射到网格上。这种复杂的行为用一组相互作用/卷积矩阵Mk，0≤k≤Kmax来表示。目的是使用卷积矩阵Mk和向量ln计算未知数a n在时间n的状态，向量l n描述在时间n在网格的未知数上发射的入射波。n，通过KmaxM k·an−k= l n.（一、一）k≥0向量a和l的维数为N，矩阵M的维数为N × N，其中N为未知数n。这是一个KmaxMk·an−k表示不使用由向量an−k（k>0）定义的未知数的预存状态的情况下，预存预处理项。这个方程表达了这样一个事实，即现在在时间n（an）发生的事情是过去n-k（an-k）发生的事情的后果6在远离k·c·Δt的元素上。这种距离关系在矩阵Mk中得到了解释。原始方程（1.1）可以重写为公式（1.2），其中左手侧是要计算的状态，右手侧是从Kmax先前时间步长和测试用例定义中已知的。对于T次迭代/时间步长n>0，我们计算an=（ M0）−1（ln−Kmaxk=1Mk·an−k）.（一、二）1.2.2相互作用/卷积矩阵矩阵Mk包含由k·c·Δt附近的距离分开的未知量之间的相互作用，并且对于比该距离更近或更远的未知量包含零它们具有以下特性：• 给定矩阵Mk的非零值的数量取决于网格的结构（未知数之间的距离）和系统的物理性质c，λ和Δt。在实际构形中，矩阵M k是稀疏的，因为远离k·c·Δt的元素只是网格的一小部分。• 当k> Kmax = 2 + Kmax/（cΔt）时，矩阵Mk为零，其中Kmax为物体的直径（Kmax=max（x，y）∈Γ×Γ（|x-y|））。在时间t发生的事件不会影响日期（t + L/c）之后的结果。• M0几乎是对角线的，因为它代表了元素与它们自己和它们的近邻的相互作用。• 非零值在矩阵中的位置由未知数的编号决定，通过适当的编号，我们希望矩阵看起来像图1.3。M0M 1M 2图1.3：Mk矩阵的期望形状...MKmax图1.4说明了这些矩阵的构造，并显示了非零值（NNZ）的位置取决于一个未知数发射的波经过另一个未知数所需的延迟两个NNZ值在一行中是连续的，例如Mk（i，j）和Mk（i，j +1），如果未知数j和j+1都在距离k处。C. Δt从i.另一方面，如果未知数i和i+1都在距离k处，则两个NNZ值在列中是连续的，例如Mk（i，j）和Mk（i+1，j）。C. Δt∑7∑·˜来自j.因此，对空间上接近的未知数进行连续编号是增加在交互矩阵中具有连续值的机会的一种方式然而，除了少数例外，不可能找到一个编号的未知数，导致完美的NNZ对角线的三维问题。A B CA B CABCA B CA B CABC(a) M0(b) M1(c) M2(d) M3图1.4：二维中三个未知数A，B，C的Mk矩阵示例.在每个时间步长处表示从每个未知量发射的波当一个波围绕一个未知数时，一个值被添加到矩阵中，该矩阵由灰色正方形表示所有k>3的矩阵Mk都是零，因为元素之间的最大距离≤3cΔt。在我们的模拟中，网格是固定的（在时间上是恒定的），这为矩阵和解决问题的方式提供了重要的属性;矩阵在每个时间步都是相同的，它们是有限数量的，并且在每个时间步都求解相同的矩阵M0。1.2.3分辨率算法解分两步计算 在第一步中，使用p和

下载后可阅读完整内容，剩余1页未读，立即下载