鹌鹑：用Python实现的轻量级不连续Galerkin求解器

软件X 17（2022）100982原始软件出版物Quail：Python中用于教学和原型设计的轻量级开源不连续Galerkin代码Eric J. Ching，BrettBornhoft，Ali Lasemi，Matthias Ihme斯坦福大学机械工程系，斯坦福，CA 94305，美国ar t i cl e i nf o文章历史记录：2021年1月25日收到2021年12月2日收到修订版，2022年保留字：间断Galerkin方法Pythona b st ra ct在本文中，我们提出了鹌鹑，一个轻量级的不连续Galerkin求解器编写的Python。本守则的目的不仅是作为这个迅速发展的领域的新来者的教学工具，而且还作为用于测试算法、物理模型和不连续Galerkin框架中的其它特征的原型平台。代码的可读性、模块化和易用性都得到了强调。目前，Quail在1D和2D非结构化网格上求解一阶和二阶偏微分方程包括各种时间步进方案、求积规则、基础类型、方程组和其他功能。代码的结构和功能，以及代表性的例子，涉及传播的二维等熵涡和二维黎曼问题的重力源项，将进行讨论。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-21-00020Code Ocean compute capsuleGNU通用公共许可证GNU General Public Licensev3.0使用git的代码版本控制系统使用Python 3.7的软件代码语言、工具和服务编译要求，操作环境依赖NumPy [1]，Matplotlib [2]，SciPy [3]如果可用，链接到开发人员文档/手册https://github.com/IhmeGroup/quail/blob/main/docs/documentation.pdf问题支持电子邮件bornhoft@stanford.edu1. 动机和意义偏微分方程的数值解与许多科学和工程领域有关，包括流体力学、固体力学、电动力学和天体物理学。各种数值方法，如有限差分和有限体积格式，可以用来离散描述这些问题的控制方程。近年来，高阶间断Galerkin（DG）方法得到了相当大的关注[4，5]。这些方法结合了经典的有限体积法和有限元法。整体解的近似使用分段不连续多项式，导致在每个元素的多个自由度。元素间的不连续性用数值通量来解释*通讯作者。电子邮件地址：bornhoft@stanford.edu（Brett Bornhoft）.https://doi.org/10.1016/j.softx.2022.100982功能协调发展的DG方案的优点包括高阶精度（通常定义为大于二阶）、理想的耗散和色散特性、几何灵活性以及适用于hp-自适应的能力，其中h是指网格，p是指精度的阶数。此外，DG方法可以在高性能计算系统上实现非常好的可扩展性和效率在空气动力学[6，7]、多相流[8，9]、等离子体物理学[10，11]、天体物理学[12，13]和固体力学[14，15]中已经证明了令人鼓舞的性能。然而，DG方案和相关的高阶方法，例如通量重建[16]和谱差分方案[17，18]，通常比低阶方案更不鲁棒并且更占用内存。此外，这些高阶方法通常需要的弯曲网格很难生成，并且可能需要非平凡的扩展来解释额外的物理学。这些因素目前阻碍了DG方案在工业应用中的广泛使用2352-7110/Elsevier B. V.这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxEric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009822��ˆ分别表示关于车身的内部和外部信息∑e+e−He=1HHH∑e+HHHH∑被认为是更复杂的学习和实施比∑矩阵MNeHH∫=∫与DG计划相关的另一个困难是，其中，Me=εφmφndε表示局部元素 质量传统的低阶方法[19]。这种进入壁垒不仅会阻碍工程师和科学家进入该领域，而且还会阻碍实施和算法的进展在方程的第二项上进行分部积分（4）、屈服胸腺发育可靠的学习资源可以帮助降低φm这个屏障。有很好的教科书可以描述这个理论，和DG方案的应用[5，20]。也不在少数H H埃莱e+e−+DG和相关高阶方法的开源代码[2126]。然而，这些通常是大型的生产级代码，F（U得双曲余切值.，n）dr，（6）阿吉耶这对新来者来说可能是压倒性的，并运行.Hethhaven和Warburton的教科书[5]包括一组用MATLAB和C++编写的著名例程，但很难以模块化的方式添加新功能在本文中，我们提出了鹌鹑，一个轻量级的DG代码设计，以解决上述问题。代码是用Python写的因为它是开源的（不像MATLAB），文档齐全，其中n是指向外的单位法向量，（·）+和（·）-eF是数值通量函数。合并等式（4）、（5）和（6）给出了以下半混凝土形式：Nb适用于面向对象编程，越来越受欢迎，e e e并且比低级语言（如C、C++和Fortran）更易于阅读和使用。鹌鹑的主要目标是作为一个教学和原型工具。为了实现这一目标，我们专注于代码的模块化、清晰度以及简单的设置和使用。Quail的简单性允许用户轻松地阅读，理解，n=1dtUn（t）Mmn埃莱--φm·F（Uh）dφφmF（U得双曲余切值.，n）dΓ并添加到代码中，而不需要大型代码库的复杂性。本文的其余部分组织如下。第2节简要概述了离散化。第3节描述了软件的结构和功能，随后是一个由二维等熵涡流传播组成的说明性示例。最后两节分别讨论这一守则的影响和结论性意见。2. 数学背景Quail求解了以下形式的线性和非线性偏微分方程组：U+F=S，（1）其中U是状态变量的矢量，F是通量，S是源项。为了简洁起见，本节假设一阶偏微分方程，尽管Quail也可以求解二阶偏微分方程。设计算域为N，不重叠的离散元素，使得 你好阿吉耶=φmS（U e）d.（七）埃莱LHS上的时间导数可以用经典的显式时间步进格式来处理。使用数值积分的积分进行评估关于离散化的更多细节可以在文档（代码元数据）中找到，以及参考。[5]的第10段。3. 软件描述Quail使用面向对象的编程范式来实现灵活性。因此，当实现一个新特性时，继承和构建现有的类结构是很简单的可读性是通过采用Python编码风格和广泛使用Python文档字符串和注释来实现的所有的数学运算都是使用流行的NumPy库[1]执行的，它通过充分利用矢量来提高性能是元素的边界e1e.近似于全局torized数组操作。 为了实现模块化，Quail被分为五个模块。解UhU可以扩展为Uh=NeUe，（2）其中Ue是局部离散解，Python包，每个包负责封装不同的DG解决方案程序的各个方面。该软件体系结构如图1所示。 1，具有表1中列出的特定功能。每个包装的简要说明如下：Ue（x，t）=eNbn=1n（t）φn（x）.（三）网格包处理计算网格及其数据结构的生成。支持高阶曲面元素。内置的例程提供生成1D和2D规则网格与四边形和un（t）是基系数的第n个向量，φn是基系数的第n个向量。基多项式基础的常见选择包括拉格朗日多项式和勒让德多项式。光滑解的空间精度的名义阶为p1，其中p是多项式近似的阶。为了求解局部离散解，我们要求Ue满足三角形元素也可以在用户指定的边界上施加周期性。该软件包还可以从开源网格软件Gmsh[27]导入网格，该软件可以为复杂几何形状生成弯曲的非结构化网格。• numerics包负责提供阿吉耶φmtUed+阿吉耶φm·F（Ue）d=阿吉耶φmS（U e）dφm.（四）数值方法这些措施包括各种选择的二次积分规则评估积分，基函数的解决方案和几何近似，时间步进计划，和限制器和人工粘性稳定。的插入等式（3）在Eq.中的第一项（4）让我们写NbNbe e e e==·˜阿吉耶公司简介∫Eric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009823求积规则和基础类型支持1D中的线段以及2D中的三角形和四边形，并且所有这些形状都有节点和层次基础φmtUhd埃莱n=1 dtUn（t）φmφnd埃莱n=1 dtUn（t）Mmn，（五）一 Gauss-lobatto 并置 方案， 在 这 解决方案分节点和求积点是相同的，包括.Eric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009824=-图1.一、 软件体系结构图描述了每个包执 行 的 任 务 以及代 码 的总体流程。物理软件包包括方程组和相应的物理模型。它包含通量、初始条件、边界条件和源项的定义。支持的方程组包括常数标量的advection，无粘Burgers方程，和可压缩的欧拉方程。我们最近还增加了处理二阶偏微分方程（如标量平流）扩散和可压缩求解程序包将存储在物理中的物理相关信息与数值算法相结合，以在每个时间步更新解决方案。 该解不仅可以用标准DG离散化计算，此外，ADERDG方案[28，29]也是一种空时预测-校正方法，允许在空间和时间上具有高阶精度，同时保持刚性源项的鲁棒处理包对解决方案数据执行任务，例如计算错误和生成可视化。包括1D线图、2D等高线图、线探测和动画。Matplotlib [2]子例程用于方便设置、实现和修改。用于处理和模拟重启的解决方案数据文件使用pickle模块保存，该模块是Python标准库的一部分。简单的子程序允许轻松地读取和写入数据。整个Python对象可以用一行代码保存，这是该模块的一个主要优点。Quail的另一个用户友好的特性是可以选择定义自定义函数，以在每个时间步执行特定于案例的数据处理这些自定义函数是在源代码之外编写的，因此不需要完全理解代码库，允许简单的实现。还包括对持续集成的内置支持和包含功能测试和单元测试的测试可以在存储库中找到一套测试用例。这些包括一维阻尼正弦波、高斯脉冲的二维平流、Sod激波管问题、移动激波、凸块上的定常流、定常无粘第4节中的示例也可提供。其他学习工具包含在Quail中。例如，线段、四边形和三角形的基函数可以被可视化，如图2所示。二、 此外，可以执行各种多项式阶数的耗散和色散分析。为了说明，图3给出了p1到p7在逆风通量下的耗散和色散关系。有关该分析的更多信息，请参见参考文献。[5]的第10段。4. 说明性实例4.1. 等熵涡传播为了说明Quail的功能，我们给出了由可压缩欧拉方程控制的传播等熵涡的数值解。这是检验高阶格式精度的一个经典例子。输入甲板依赖于Python字典来规定求解器和物理参数，允许轻松的模拟设置。这些输入字典被组织如下：数学，网格，物理，初始条件，精确解，BoundaryConditions、SourceTerms、Output和Restart。位于src/defaultparams.py中的默认参数只允许用户修改每个人的必要参数案子由于Python脚本被用作输入文件，因此任何Python功能都可以直接在输入组中使用。这种情况下的输入文件如下所示#用于设置#propagatingisentropic的输入卡组vortexTimeStepping={“FinalTime“：1. 0、“CFL“：0.1，“TimeStepperr“：“LSRK4“，}Numerics={“SolutionOrder“：3，“SolutionBasis”：“LagrangeTri“，}Mesh={“ElemmmtShape“：“Triangle“，“NumElemsX“：16，“NumElemsY“：16，“xmin“：-5. 、···Eric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009825[−][−] ××=我的天=== − −=图二、（ a）线段和（b）四边形的样本基函数。邮政局 get_error（mesh，physicss，solver，“Density“）#Plotdensity 我的天啊plot。 prepare_plot（inewidth =0. 第一章plot。plot_soluton（mesh，ph ysics，solver，“Densi ty“，lege nd_label=“DG“，inclde_mesh=True，reg ular_2D= True）plot. save_figure（文件e_name=）plot。show_plot（）初始条件基于Yu等人的工作。[31]，其中速度和温度扰动叠加在均匀流上。域是大小为5的正方形，五五精确解对应于伊森的平流，恒速的热带涡旋。图图4示出了在两个不同网格上在一秒的最终时间此脚本读取pickle数据文件，从求解器中解包相关对象，计算密度的L2误差，并绘制密度等值线。这些数据可以用来绘制收敛图，如图11所示。4（d），这表明获得了预期的收敛速度。有关如何构造输入文件和后处理脚本的其他信息，包括代表性示例，可以在Quail存储库（代码元数据）中找到。4.2. 带重力源项的二维Riemann问题Quail的面向对象框架使得添加源项和数值算法变得简单。在这里，我们说明了使用这些功能，通过模拟一个2D黎曼问题的重力源项由可压缩欧拉方程。如参考文献[32]中所做的那样，在具有以下初始条件的2 2域上建立该情况：1个三角形网格（图4（a）），其中q是指的几何近似，和一个弯曲的q2四边形网格（图。4（c））。模拟完成后，求解器（ρ，u，v，P）（7，−1，0，0. 2），x≤1，（7，1，0，0. 2），x> 1，（八）（可选）在工作目录中搜索后处理脚本。使用第3节中详细介绍的工具，用户可以创建等高线和线图，如图1所示。 四、本例的后处理脚本如下：我很抱歉。PostasPost我很抱歉。 plotasplot我很抱歉。读取已写入的数据文件#读数据文件fname=“Data_final. pkl“solver=read写了一个文件。read_data_file（fname）#Unpackmesh=solver.物理学=求解器。物理#ComputeL2errorr其中ρ是密度，u和v是x和y分量P是压强。 源项表示为S[0，0，ρg，ρv g] T，with g1. 在这些条件下，压力和密度都可以接近非物理负值.为了提高鲁棒性，我们采用Zhang和Shu [33]的保正限制器。多项式的阶数是p1。图 5，我们给出了在t处的Riemann问题的数值结果0的情况。6s，这是可比的那些在文学[32]。5. 影响虽然DG方法越来越受欢迎，但由于它们通常更复杂，因此“xmax“：5.，“ymin“：-5. 、“ymax“：5.，}Physics={“类型“：“Euler“，“ConvFluxNumerical“：“Roe“，“GasConstant“：1. 、}InitialCondition={“Function“：“IsentropicVortexx“，}ExactSolution=InialC onion。publicintfindDuplicate（）{“BCType“：“StateAll“，“Function“：“IsntropicVortex“，}在s={“x1“：d，“x2“：d，“y1“：d，“y2“：d，}Eric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009826图3.第三章。（ a）耗散和（b）p = 1至p = 7的 色 散 关 系 ，有逆风通量。图四、等熵涡传播的数值解。（a）密度等值线图a tt=1。1024个三角形元素上的0。（b）y=1处的线切片比较初始、精确和数值解。 (c)收敛速度的变化基于L2密度误差的不同多项式阶数的Ne期望的收敛速率用实线表示Eric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009827==表1Quail中实现的当前功能集基础与几何形状段三角形四边形节点基础拉格朗日拉格朗日拉格朗日模态基础LegendreLegendreH1层次[34]求积规则Gauss–Legendre,Gauss–Legendre,Gauss–Legendre,稳定时间步进器[33]第三十三话[36]第三十六话人造粘度LSRK 4 [37][38]第三十八话DG ADER [39][40]第四十话：一个人的世界[41]第四十一话标量图五. 带重力源项的二维黎曼问题的数值解。(a)t时的密度等值线0的情况。25，600个四边形元素上的6 s。(b)线在y 1处切片。7875比较6400个元素情况（蓝色符号）与25，600个元素情况（黑色实线）的数值解。(For对本图图例中所指颜色的解释，读者可参考本文的网络版比传统的低阶数值格式，如有限体积法和有限差分法，有相对较少的学习资源。为了帮助解决这个问题，鹌鹑旨在促进学生，科学家和渴望进入该领域的工程师的学习。在文献中不容易讨论的关键实现细节通过这个轻量级的用户友好的与大型代码库的复杂性相比，Quail的简单性使其更容易被那些希望亲自访问DG计划的新手所消化。这种简单性，以及模块化和Python作为首选语言，使Quail也有利于快速原型设计。物理模型可以很容易地纳入评估性能的DG计划在新的背景下。此外，新的功能和方法可以在添加到生产代码之前快速实现，测试和应用于模型问题。这可以加速算法开发的过程。矢量化NumPy操作的广泛使用[1]使得能够模拟具有相关物理尺度和分辨率的各种1D和2D配置。存储库中提供的视频教程链接（代码元数据）进一步说明了如何使用Quail并以简单的方式添加某些功能，例如边界条件。6. 结论Quail是一个用Python编写的轻量级不连续Galerkin代码。它是专为教学和原型设计的，没有生产代码的笨重复杂性。代码清晰、模块化和易用性是主要的关注点。目前，Quail解决了1D和2D一阶和二阶偏微分方程。软件架构和功能进行了讨论。等熵涡的传播和一个二维黎曼问题的重力源项由欧拉方程作为说明性的例子。Quail可以通过降低新来者的进入门槛和加速算法开发来影响社区。未来的工作将需要纳入新的功能，如额外的限制器，方程组和数值计划。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作Eric J.Ching，Brett Bornhoft，Ali Lasemi等.软件X 17（2022）1009828致谢来自NASA空间技术研究资助计划，国家科学基金会（DGE-1656518和1909379）和国防部SMART奖学金的早期职业教师资助（NNX 15 AU58 G）表示感谢。引用[1] Harris CR，Millman KJ，van der Walt SJ，Gommers R，Virtanen P，Cour-napeau D，Wieser E，Taylor J，Berg S，Smith NJ，KernR，皮库斯M、Hoyer S，van Kerkwijk MH，Brett M，Halflian A，del Río JF，Wiebe M，Peterson P ， Gérard-Marchant P ， Sheppard K ， Reddy T ， Weckesser W ，Abbasi H，Gohlke C，Oliphant TE.使用NumPy进行数组编程。Nature2020;585（7825）：357-62.[2] 亨特JD。Matplotlib：2D图形环境。计算机科学与工程2007;9（3）：90-5.[3] Virtanen P，Gommers R，Oliphant TE，Haberland M，Reddy T，CournapeauD，Burovski E，Peterson P，Weckesser W，Bright J，van der Walt SJ，BrettM，Wilson J，Millman KJ，Mayorov N，Nelson ARJ，Jones E，Kern R，Larson E，Carey CJ，Polat I，Feng Y，Moore EW，VanderPlas J，LaxaldeD ， Perktold J ， Cimrman R ， Henriksen I ， Quintero EA ， Harris CR ，Archibald AM ， Ribeiro AH ， Pedregosa F ， van Mulbregt P ， SciPy 10Contributors.SciPy 1.0 ： Python 中 科 学 计 算 的 基 本 算 法 。 NatureMethods2020;17：261-72.[4] Cockburn B， Karniadakis GE ，Shu CW.间 断 Galerkin 方 法的 发 展。 In ：Cockburn B ， Karniadakis GE ， Shu CW ， editors. 间 断 Galerkin 方 法Springer，Berlin，Heidelberg; 2000，p. 三比五十[5] Hesthaven JS ， Warburton T. 节 块 间 断 伽 辽 金 法 ： 算 法 、 分 析 与 应 用 。Springer; 2007.[6] 刘 晓 波 ， 李 晓 波 ， 刘 晓 波 .计 算 空 气 动 力 学 的 间 断 伽 辽 金 法 - - 用 D L RP A D G E 代 码 进 行 三 维 自 适 应 流 模 拟 。 Aerosp Sci Technol 2010;14（7）：512[7] Brazell MJ，Mavriplis DJ，Yang Z. 使用有限体积和间断Galerkin解算器的网格分辨翼型模拟。AIAA J2016;54（9）：2659-70。[8] Zwick D ， Balachandar S. 快 速 减 压 多 相 激 波 管 的 动 力 学 。 J FluidMech2019;880：441-77.[9] 张文辉，张文辉，张文辉.用双向耦合欧拉-拉格朗日方法在任意曲面单元上用间断伽辽金格式模拟多相流。第405章.[10]Jacobs G，Hesthaven JS.非结构网格上的高阶节块间断Galerkin质点网格法。JComput Phys2006;214：96-121.[11]张文辉，张文辉，张文辉，张文辉，张文辉. 基于高阶杂交间断Galerkin谱元法的非结构曲面网格的粒子网格求解器。计算方法应用机械工程2019;349：149-66。[12]Chu R，Endeve E，Hauck CD，Mezzacappa A.费米子输运两矩模型的保可实现性DG-IMEX方法。J ComputPhys2019;389：62-93.[13]Teukolsky SA.相对论天体物理学的间断伽辽金方法的公式化。J ComputPhys2016;312：333-56.[14]作者：J. H.非线性弹性力学的可杂交不连续Galerkin公式。计算方法应用机械工程2015;283：303-29.[15]Nguyen VP.不连续Galerkin/非本征内聚区模型：计算断裂力学中的实施注意事项和应用。工程机械2014;128：37-68。[16]惠恩海峡含间断Galerkin方法的高阶格式的通量重构方法。第18届AIAA计算流体动力学会议。AIAA 2007-4079;2007年。[17]刘毅 ，范诺 科，王 正杰 .非结 构网格 的谱差分 法I：基 本公式 。J ComputPhys2006;216（2）：780-801.[18]王志杰，刘燕，梅G，詹姆森A.非结构网格的谱差分法II：扩展到欧拉方程。JSci Comput2007;32（1）：45-71.[19]Wang ZJ，Fidkowski K，Abgrall R，Bassi F，Caraeni D，Cary A，DeconinckH，Hartmann R，Hillewaert K，Huynh H，Kroll N，May G，Wavesson P-O，vanLeer B，Visbal M.高阶CFD方法：现状与展望。Int J Numer MethodsFluids2013;72：811[20]里维耶尔湾解椭圆和抛物方程的间断伽辽金方法：理论与实现。SIAM; 2008年。[21]坎特韦尔CD，莫克西D，科默福德A，博利斯A，罗科G，门加尔多G，DeGrazia D ， Yakovlev S ， Lombard J-E ， Ekelschot D ， Jordi B ， Xu H ，Mohamied Y，Eskilsson C，Nelson B，Vos P，Biotto C，Kirby RM，SherwinSJ. Nektar ++：一个开源的spectral/hp元素框架。Comput PhysComm2015;192：205-19.[22][10]杨文，李文，李文，李文. 非定常问题的显式间断Galerkin方法。Comput&流体2012;61：86[23]Alnæs M，Blechta J，Hake J，Johansson A，Kehlet B，Logg A，RichardsonC ，Ring J ， Rognes ME ， Wells GN. FeNiCS 项 目 版 本 1.5 。 ArchNumerSoftw2015;3（100）.[24]Witherden FD，Farrington AM，Vincent PE. PyFR：一个开源框架，用于使用 通 量 重 建 方 法 解 决 流 架 构 上 的 对 流 扩 散 类 型 问 题 。 Comput PhysComm2014;185（11）：3028-40.[25][10]杨文辉，陈文辉. The Deal.II有限元库：设计、功能和见解。计算数学应用2021;81：407[26]Klöckner A，Warburton T，Hesthaven JS.在非结构几何上求解波动方程.在：GPU计算宝石翡翠版. Elsevier;2012，p. 225比42[27]Geuzaine C，Remacle JF. Gmsh：一个内置预处理和后处理设施的三维有限元网格生成器。Int J Numer Methods Eng2009;79（11）：1309-31.[28]Dumbser M，Zanotti O.非结构网格上电阻性相对论MHD方程的极高阶P N PM格式J Comput Phys2009;228：6991-7006.[29]作者：Bornhoft BJ，Ching EJ，Ihme M.用间断伽辽金法解反应流的时间积分考虑。在：AIAAScitech2021论坛。AIAA 2021-0745; 2021.[30]王丙重构间断Galerkin方法求解任意拉格朗日和欧拉形式的可压缩Navier-Stokes方程（博士）论文），北卡罗来纳州立大学; 2017年。[31]余明，王正，刘勇.关于间断伽辽金法的精度和效率，谱差分法和重建法的修正过程。J Comput Phys2014;259：70-95.[32]张 X ， 舒 C-W. 含 源 项 可 压 缩 Euler 方 程 的 保 正 高 阶 间 断 Galerkin 格 式 JComputPhys2011;230（4）：1238-48.[33]张X，舒C-W.矩形网格上可压缩Euler方程的保正高阶间断Galerkin格式。JComput Phys2010;229（23）：8918[34]作者：J. M.高阶有限元素法。CRC;2003，p. 55比60[35]杜纳万特湾三角形的高效率对称高斯求积规则。Int J Numer Methods Eng1985;21：1129[36]Zhong X，Shu C-W. Runge-Kutta间断Galerkin方法的一个简单的加权本质无振荡限制器。J Comput Phys2013;232（1）：397-415.[37]Carpenter MH，Kennedy C.四阶2N存储Runge-Kutta格式NASA兰利研究中心，NASA报告TM 109112，里士满，弗吉尼亚州，美国;1994年。[38]Spiteri R，Ruuth S.一类新的最优高阶强稳定性保持时间离散方法。SIAM JNumer Anal2002;40（2）：469-91.[39]Dumbser M，Enaux C，Toro E. 刚性双曲平衡律的高精度有限体积格式。JComput Phys2008;227：3971-4001.[40]斯特朗湾差分格式的构造与比较。SIAMJ Numer Anal1968;5（3）.[41]吴H，马P，Ihme M. 模拟具有刚性化学反应的湍流反应流的有效时间步进技术。Comput Phys Comm2019;243：81-96.