ExWave：高性能声波方程间断Galerkin求解器

SoftwareX 9（2019）49ExWave：一种高性能的声波方程间断Galerkin求解器S. Schoeder，W.A.沃尔，M。克龙比希勒德国慕尼黑工业大学计算力学研究所ar t i cl e i nf o文章历史记录：收到2018年2019年1月2日收到修订版，2019年MSC：65M6065Y2068N19保留字：无矩阵方法间断Galerkin方法声波方程a b st ra ct本文提出了一种高性能的间断Galerkin离散方法，并采用显式对于ADER，提供全局和局部时间步进变量该实现是基于无矩阵框架的deal。II有限元库提供了高效的评估程序的四边形和六面体。该实现是通用的，其适用性证明了学术的例子，以及现实世界的问题，如城市声学。我们提出的物理和数值问题的描述，一般的代码结构，和设计原则。©2019作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0此代码版本使用的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2018_172法律代码许可证LGPL v.2.1使用Git的代码版本控制系统使用C++、MPI的编译要求，操作环境&依赖cmake，ctest，deal.II version 9.0/9.1，p4est如果有开发者文档/手册github的链接。com/kronbichler/exwave问题支持电子邮件kronbichler@lnm.mw.tum.de1. 动机和意义声波方程有着广泛的应用.除此之外，它还可以预测建筑项目的室内声学效果，或优化城市规划和城市设计的噪音。存在无数的数值工具来求解声波方程，每一种都有其优点和缺点。最具挑战性的是高频波的精确模拟，因为它们需要高的空间和时间分辨率。一组解决方案技术规避了由高频引起的困难，并适用于室内声学是基于射线跟踪的几何方法，然而，这对于低频和衍射不够准确[1]。另一种方法是假设声波在房间中扩散传播，*通讯作者。电子邮件地址： schoeder@lnm.mw.tum.de（美国）Schoeder）。网址：http://www.lnm.mw.tum.de（S. Schoeder）。https://doi.org/10.1016/j.softx.2019.01.001足够数量的反射和扩散方程模型被求解（例如，参见[2]）。在很宽的频率范围内的声学预测需要用具有足够精细离散化的纯数值方案或部分依赖于波动方程的基本解的半解析方案来求解声波方程。简单的时域有限差分（FDTD）方法被广泛使用，但由于有限的计算资源，仅适用于较低的频率[3]。最近，已经提出了自适应矩形分解， 区域分解技术依赖于矩形波方程的解析解和附加的界面处理[4，5]。它们目前仅限于均匀声速分布。除了有限差分格式外，常规有限元、混合元、谱元、间断Galerkin（DG）方法和有限体积方法也分别应用于声波方程，见[62352-7110/©2019作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx50S. Schoeder，W.A.Wall和M.Kronbichler/SoftwareX 9（2019）49[+[客户端]S∑∑Pti（）+=[]个=阿[客户端]∈（[]表1实施时间积分方案。所示订单指的是时间积分方案。对于ADER，也可以达到k+2阶，如果向量V、P中总结的未知数。对象的更新规则s-阶段Runge-Kutta格式如[12]中的重建名称描述订单阶段Explicit EulerExplicit Euler1 1Vti1Pti+1VtiPti公司简介∑j=1 bjKj，其中（3）clRK4经典RK 4 4LSRK45R2低存储RK，两个寄存器[17] 4 5LSRK45R3低存储RK，三个寄存器[17] 4 5LSRK33R2低存储RK，两个寄存器[17] 3 3Kj= −Q−1KVtiPti公司简介j−1l=1ajlKl），LSRK59R2低存储RK，两个寄存器[17] 5 9SSPRKStrong-stability preserving RK [18] 4 8ADERADER如（4）或[12]k+ 1中所示。ADERLTSADER LTS，如[12]k+ 1-ADERADCONFULL伴随一致ADER [12]k+1-以及来自相应方案的Butcher表的系数a，j， b，j。其中，λt是时间步长，矩阵Q和K是由空间离散化产生的质量和刚度矩阵，如[12]所示对于ADER时间离散化，更新规则为这项工作的基础是关于DG的出版物[11，12][Vti+1][Vti]-Q−1KQ−1K+1tj+1（−1）j声波方程的方法。在这篇文章中，我们提出了一个高性能的高阶DG求解器的声波方程，是一般的几何形状，频率，ti+1我不是×∫Ωj=0 （j+1）！NTSjNd[Vti]，（4）范围和声速分布，并与显式Runge-Kutta以及任意导数（ADER）时间步进方案相结合我们的代码支持空间中的自适应网格细化以及时间中的本地时间步进，以将计算工作集中在最有趣的物理现象发生的区域。该实现基于deal。II有限元库[14]，为四边形和六面体[15，16]提供有效的无矩阵评估例程该算法在吞吐量、每秒处理的自由度和可扩展性方面在现代硬件上提供了非常高的性能[13]。ExWave对科学计算社区很有价值，因为它平衡了高性能和代码优化与现实世界问题和可用硬件的适用性。对于计算声学社区，该代码允许接近当前未探索的场景，因为其在几何形状，网格生成或自适应性方面具有很高的灵活性（如[12]），但也可以在合理的时间内以适度的计算资源提供可预测精度的解决方案。理论背景以压力p和粒子速度v表示的写为一阶系统的声波方程给出为：其中k是所用形状函数的多项式次数，N是保持形状函数的矩阵，算子S包括表示波动方程的空间导数，详见[12]。方法（3）和（4）都给出了在压力p和速度v下k1阶的最优收敛。通过重构和后处理，甚至可以得到p中k+2阶的超收敛结果数值算例虽然实现是相当普遍的（参见例如，第3.3节中的示例），我们将结合一个学术示例对其进行解释，该示例的参数文件在当前代码版本中提供例子是振动膜问题的解析解是可用的，并实施，因此允许简单的收敛性测试。在d维立方体上，每个坐标有m个振型的振动膜服从精确解Dp=cosMdπt·sin（mπxe）（5）e=1齐次Dirichlet边界条件对更一般的设置的修改可以通过如3.3节所示的代码的扩展来容易地引入。吉夫公司简介1ρρ=0，（1）参数文件允许调整空间维度（d2，3）、空间离散化（即，网格、网格细化、网格变换）。此外，关于速度的参数np+c2ρε·v=0，（2）声速c和质量密度ρ。这个方程在d维空间域上成立，研发和在时间间隔0，T与最终时间T。第一个方程表示动量守恒，而第二个方程强制执行质量守恒。声波方程伴随有关于压力和速度的初始条件以及边界条件。目前，在ExWave中实现了压力的Dirichlet边界条件、速度法向分量的齐次Neumann边界条件如[11]所述，使用来自混合DG方法的精确通量，通过DG方法进行空间离散化。时间离散化由显式Runge可以调整Ral离散化。对于时间积分，可以在表1中列出的方案之间进行选择。此外，用户输入最终时间、输出时间步长和导出时间步长方程中的参数m（5）是由membrane_modes指定。 除此之外，可以使用use_ader_post设 置 ADER 和 ADER LTS 特 定 参 数 ， 启 用 重 建 和spectral_evaluation，参考[ 13 ]中介绍的快速评估。2. 软件描述所提供的代码基于deal.II有限元库的9.1-pre1版本[14]。为了确保代码的功能性和正确性，并简化进一步的开发，使用ctest测试框架结合了检查几个设置的收敛顺序的单元测试。Kutta格式，如[11]中所示，或通过ADER时间步进，在[12，13]。 该解决方案表示为2018年11月28日102d808之后的时间依赖1 Any commit。PP]=S. Schoeder，W.A.Wall和M.Kronbichler/SoftwareX 9（2019）4951×+=·+=+=-Fig. 1. 主类WaveEquationProblem包含一个用于计算空间离散化WaveEquationOperationBase的运算符和一个来自ExplicitIntegrator的时间积分器。 阴影区域表示如何组合空间和时间运算符。ExWave的主要类是WaveEquationProblem。 它的方法run（）执行时间循环。主要组件是从ExplicitIntegrator导出的时间积分器和从WaveEquationOperationBase导出的空间运算器，如图所示。1.一、时间积分器执行向量更新并调用空间运算符的应用。对于ADER，空间和时间评估紧密相连，整个评估在WaveEquationOpeationADER中进行。LTS需要COM-plex更新调用，该调用由WaveEquationOperationADERLTS依次调用的DataManager处理。类WaveEquationOperation是模板化的，形状函数的阶数d和多项式阶数k。它严重依赖deal.II库的MatrixFree 类 ， 并 使 用 优 化 的 评 估 例 程 FEEvaluation 和FEFaceEvaluation。与经典的基于矩阵的计算相比，无矩阵运算符计算允许更高的性能。方案，这是由于更高的算术密度。此外，使用依赖于利用形状函数的张量积结构的和因式分解的快速积分技术，并且启用显式跨单元向量化。关于无矩阵方法的详细信息，可以在声波方程DG的上下文中找到和因子分解技术和性能[15]，[16]，[17]，[18]，[19]。WaveEquationOperationBase可以通过typedefvalue_type在单精度和双精度之间灵活切换。在文件time_integrators. h中，不仅时间积分器而且实现了优化的矢量更新器向量更新器RKVectorUpdater将多个向量更新合并到条目上的单个循环中，从而将每个阶段所需的向量读取数量从五个减少到两个。这种贡献，灰允许一个增加性能的因素1。7对于现代处理器上的LSRK45R2，其中性能通常是内存带宽有限的[13]。我们要明确提到的代码的一个方面是，它允许运行一个迭代CFL稳定性分析的基础上的稳定性标准的L2压力误差，找到一个紧配合的临界柯朗数为一定的问题配置。3. 说明性实例在这一章中，给出了一个学术示例，然后是两个性能测试，最后是一个现实世界问题的代表3.1. 收敛性测试通过收敛性测试证明了方法和代码的正确性为了运行基于振动膜的收敛测试，必须指定一组基本的输入参数，然后运行几个模拟，改变参数n_refinements。输出最终时刻的L2压力误差。示出了两种不同示例性配置的结果。第一种设置是二维几何形状，并且时间离散化依赖于ADERLTS。此外，设置参数use_ader_post=true以获得超收敛结果[12]。第二种设置是三维和时间离散化是基于LSRK33R2。对于这两种设置，解析解中的膜模式的数量被设置为形状函数的多项式次数，并且cfl_number=0.1，其中n_initial_intervals= 5。所有其他参数的设置如默认参数文件。图 2总结了测试多项式次数k的结果一二三四五六k的期望收敛阶1表示压力p，k2表示后处理压力p获得所有测试的多项式次数与ADER LTS。对于LSRK33R2，预期的收敛速度只获得粗离散。对于精细离散化，时间误差占主导地位，因为LSRK33R2只有三阶3.2. 绩效评价我们比较了ADER和LSRK 45 R2在2.6 GHz双插槽Intel XeonBroadwell E5- 2690 v4机器的28个核心上的吞吐量，使用g++编译器6.2版在优化级别-march=haswell-03-funroll-loops编译。设置如第3.1节所述，具有三维几何结构，其中对于k1， 2，3，803个单元和对于k403个单元进行网格化4，. . . ，12. 图 3绘制结果。吞吐量可达6. 33 108度的自由度每秒ADER与k3达到，这是3.9倍，高于相应的吞吐量LSRK45R2。对于更高的多项式次数，ADER的优势缓慢下降，这是由于ADER是k1阶的，对于更高阶有额外的计算，而LSRK45R2对于所有多项式次数都是4阶的为了证明并行能力，我们进行了一个强大的缩放实验与h-适应性的基础上的三维52S. Schoeder，W.A.Wall和M.Kronbichler/SoftwareX 9（2019）49=-××·=·图二、ADERLTS二维收敛性研究和LSRK33R2三维收敛性研究。图三. ADER和LSRK45R2的d3和k1的测试结果表明，. .，在一个28核的节点上有12个。如[13]第4.7节中的设置，结合了阻抗失配为4的材料不均匀性和利用高阶映射捕获两种材料之间的弯曲界面的非结构化自适应更新每500个时间步长进行图4显示了在SuperMUC Phase 2系统上获得的结果，该系统使用112至7168个Intel Haswell E5-2697 v3内核，运行频率为2.1 GHz。计算时间几乎理想地缩放，而如[19-21 ]中所解释的那样，可以观察到自适应的放缓，并且由于deal.II中的h-自适应性的弱标度是更好的，因为自适应性的主要成本是由于交易的成本。II的通用传输程序之间的网格水平，几何评价在幽灵细胞，和成本的数据传输，而重新分区。 写入输出在少量内核时扩展良好，但在1792个内核后会下降。3.3. 城市声学该示例检查了代表室外声学的村庄中的声音传播，这与城市规划和城市设计相关[22]或在犯罪控制背景下对枪声定位有用[23]。所考虑的几何形状是基于[24，25]中提出的人工村庄，其中FDTD方法和自适应矩形分解方法是图四、完整模拟的强大缩放功能，具有43，000个时间步长和高达1.8亿个空间未知数，ADER为d=3和k= 5。图五. 几何训练村。用于求解声波方程。该几何体由15个不同高度的建筑物 图第五幅描绘了建筑物的几何形状。计算域是一个长方体的大小175 140 14与建筑物切断。对应于墙壁，屋顶和地面的表面被假定为完全反射，而一阶吸收边界条件被应用于所有其他边界。声源位于与来自[25]的SP1相对应的（62， 104， 1）为了读取外部生成的网格并设置相应的边界条件和初始压力场，对input_parameters. h是必要的，如在分支城市声学的提供的软件。 在[25]中，在由1. 1 × 10 - 7个网格点和1× 10 -3的时间步长。85 10−4具有2000个时间步长，其对应于450 Hz的仿真频率，在单核CPU机器上花费20分钟在这里，我们跑S. Schoeder，W.A.Wall和M.Kronbichler/SoftwareX 9（2019）4953··=·=见图6。 z=1时xy平面上和x=62时yz平面上训练村中的压力快照。对平均网格间距1.2，k 3的离散化进行模拟，得到1。15 10 7个网格点，带网格点二、4510-5和19704个时间步长达到相同的最终时间。图6绘制了几个压力快照，以给出声音传播模式的印象在双插槽Intel Xeon E5- 2690 v4 Broadwell 2.6 GHz系统的28个内核上进行仿真需要2. 5 103秒，对应于2000个时间步长的76 CPU分钟因此，与自适应矩形分解相比，具有相同数量的网格点和时间步长的相同几何形状上的计算仅慢3.8倍。考虑到自适应矩形分解依赖于在分解矩形中使用离散余弦变换的半解析方法，并且考虑到它们的时间积分不是五阶而是两步类型，这是非常好的未来的工作应解决的计算性能，考虑精度和时间稳定性的比较。4. 影响所提出的算法已在多篇研究论文中使用[11这部法律影响了两个社区。首先，对于科学计算社区，该代码是一个非常有效的实现无矩阵运算符评估的代表，其吞吐量接近更简单的有限差分方法，尽管支持复杂的网格，包括局部网格细化和弯曲网格。该代码允许运行结论性的性能测试，如[13]和3.2节所示，因此使关注高性能计算的研究人员这一新概念可以应用于声波方程以外的各种问题，如电磁学或地震学中的其他波此外，如[26]中的结果所示此外，它是引入先进功能的理想测试平台，例如用于在未拟合网格上表示复杂几何形状的切割有限元技术，例如正在进行的工作[27]。其次，这一代码将对计算产生影响，声学.历史上占主导地位的有限差分和半解析方法，我们现在有一个代码在手，这不仅是高性能，但也适用于实际相关的问题。与有限差分求解器相比，它提供了更大的灵活性，例如在网格和自适应性方面。5. 结论提出了基于高阶DG空间离散并结合显式Runge-Kutta和ADER求解声波方程的高性能程序ExWave全局和局部时间步进依赖于deal.II有限元素库，或者更准确地说，矩阵自由框架提供快速求积与和因子分解。此外，ExWave还允许在空间和时间收敛方面验证实现，在CFL稳定性分析方面确定时间稳定性限制，以及测量和比较不同离散化的计算性能。城市声学模拟，如第3.3节所示，使城市规划或街道峡谷设计最大限度地减少噪音暴露。一个潜在的扩展是实现室内和城市声学中常用的各种边界条件未来的研究应该比较这种求解器的自适应矩形分解方法，不仅在计算时间，而且精度。致谢这项工作得到了德国研究基金会（DFG）和慕尼黑工业大学开放获取出版计划的资助。此外，作者感谢德国研究基金会（DFG）通过优先计划“Ex尺度计算软件”（SPPEXA）中的“Ex尺度高阶不连续Galerkin”（ExaDG）项目提供的支持，赠款协议编号为KR 4661/2-1和WA 1521/18-1。引用[1]作者：Siltanen S，Lokki T，Savioja L.是射线还是波？了解计算室内声学建模技术的优点和缺点。在：会议记录的国际会议室声学。Melbourne，Australia;2010，p. 1比6[2]作者：John J，John J.一个扩散方程模型的局限性用 于均 匀 尺寸 房 间 的声 学 预测 J Acoust SocAm 2010;128（ 4） ：1586-9.http://dx.doi.org/10.1121/1.3479756网站。[3]波特杜伦湾低频房间声学问题的时域差分模拟。J Acoust Soc Am 1995;98（6）：3302-8. http://dx.doi的网站。org/10.1121/1.413817。[4]杨文辉，陈文辉，陈文辉.高效准确的声音传播使用自适应矩形分解。IEEE Trans Vis Comput Graphics 2009;15（5）：789-801. http://dx.doi.org/10.1109/TVCG.2009.28网站。[5]Morales N，Mehra R，Manocha D.基于矩形分解的分布式存储结构并行时域波形模拟器。ApplAcoust2015;97：104-14.http://dx.doi.org/10.1016/j.apacoust.2015.03的网站。017号[6]作者：J. 声学的有限元近似波动方程：误差控制和网格自适应。东西方数学杂志1999;7（4）：232-82.[7] Cohen G， Faugelix S.瞬 态波 动方 程 的质 量 集中 混合 有 限元 J ComputAcoust2000;8（1）：171-88.http://dx.doi.org/10.1142/S0218396X0000011X.[8]Komatitsch D，Tromp J.谱元法在三维空间中的应用三维地震波传播Geophys J Int 1999;139：806网址：//dx.doi.org/10.1046/j.1365-246x.1999.00967.x网站。[9]Dumbser M，Käser M.非结构网格上弹性波的任意高阶间断Galerkin方法-II 。 三 维 各 向 同 性 情 况 。 GeophysJ Int 2006;167 ： 319-36.http://dx.doi.org/10.1111/j的网站。1365-246X.2006.03120.x。[10]莱维克 双曲问题的有限体积法。剑桥应 用 数 学 教 科 书 ， 剑 桥 ; 2002 年 ， http://dx.doi.org/10.1017/CBO9780511791253。54S. Schoeder，W.A.Wall和M.Kronbichler/SoftwareX 9（2019）49[11] [10] Kronbichler M，Schoeder S，Müller C，Wall W.声波方程隐式和显式杂交间断Galerkin方法的比较Internat J Numer Methods Engrg 2016;106（9）：712-39. 得双曲正切值. doi.org/10.1002/nme.5137网站。[12]张文辉，张文辉，张文辉.声波方程的任意高阶显式杂交间断Galerkin方法 J SciComput2018;76 （ 2 ） ： 1-38. http://dx.doi.org/10.1007/s10915-018-0649-2网站。[13]杨文龙，李晓梅，李晓梅.波动问题高阶间断Galerkin格式的有效显式时间步长。SIAMJSciComput2018;40（6）：C803-26.http://dx.doi.org/10.1137/18M1185399网站。[14]Alzetta G， Arndt D，Bangerth W，Boddu V，Brands B，Davydov D，Gassmoeller R，Heister T，Heltai L，Kormann K，Kronbichler M，MaierM，Pelteret J-P，Turcksin B，Wells D.协议。II库，9.0版。J Numer Math2018;26（4）：173-83. http://dx.doi.org/10.1515/jnma-2018-0054网站。[15]放大图片作者：Kronbichler M，Kormann K.基于并行单元的有限元算子应用的 通 用 接 口 。 Comput Fluids 2012;63 ： 135-47. http://dx.doi 的 网 站 。org/10.1016/j.compfluid.2012.04.012。[16]放大图片作者：Kronbichler M，Kormann K.不连续Galerkin有限元算子的快速无矩阵评估，arXiv预印本（2017）http://arxiv.org/abs/1711.03590v1。[17]Kennedy CA，Carpenter MH，Lewis RM.可压缩纳维尔-斯托克斯方程的低存 储 显 式 龙 格 - 库 塔 格 式 。 应 用 数 学 2000;35 ： 177-219. http ：//dx.doi.org/10.1016/S0168-9274（99）00141-5.[18]Kubatko EJ，Yeager BA，Ketcheson DI.间断Galerkin方法的最优强稳定性保持 Runge-Kutta时间离 散。JSciComput2014;60：313-44.http://dx.doi.org/10.1007/s10915-013-9796-7网站。[19]放大图片作者：J. p4 est：paral的可扩展算法-八叉树森林的自适应网格加密SIAMJSciComput2011;33（3）：1103-33.http://dx.doi.org/10.1137/100791634网站。[20]Burstedde C，Holke J.非一致自适应网格的四面体空间填充曲线，arXiv预印本（2017）https://arxiv.org/abs/1509.04627v2。[21]杨伟华，王伟华，王伟华.大规模并行通用有限元程序的算法和数据结构。ACMTransMathSoftware2011;38（2）.http://dx.doi.org/10.1145/2049673.2049678网站。[22]康 ·J城 市 声 环 境 。 2 Park Square ， Milton Park ， Abingdon ， OXonPX 144RN：Taylor and Francis;2007.[23][10]李文，李文.城市地形中的多个同时声源定位。第四届传感器网络信息处理国际研讨会。2005年，第491-6页。http://dx.doi.org/10.1109/IPSN.2005.1440982网站。[24]作者：Albert DG，Liu L.城市环境中建筑物对声脉冲传播的影响。J AcoustSoc Am 2010;127（3）：1335-46. http://dx.doi.org/10.1121/1.3277245.[25]Morales N，Mehra R，Manocha D.城市环境中声脉冲传播的三维数值模拟。JAcoust Soc Am 2014;135（6）. http://dx.doi.org/10.1121/1.4874495网站。[26]杨伟，王伟.无矩阵高阶不连续Galerkin可压缩Navier-Stokes解算器：湍流不可压缩流的可压缩和不可压缩公式的性能比较。Internat J Numer MethodsFluids 2019;89（3）：71-102. http://dx.doi.org//fld.4683.[27]Schoeder S，Sticko S，Kreiss G，Kronbichler M.声学问题的局部时间推进高阶截断间断伽辽金方法。2018年，出版发行。