面向目标的自适应时空有限元程序的高效数据结构与算法

软件X 10（2019）100239原始软件出版物面向目标的时空有限元程序自适应性的高效可扩展数据结构和算法Uwe Köcher，Marius Paul Bruchhäuser，Markus Bause赫尔穆特-施密特大学，联邦武装部队大学汉堡，数值数学，Holstenhofweg 85，22043汉堡，德国ar t i cl e i nf o文章历史记录：2018年12月20日收到收到修订版2019年4月17日接受2019年4月17日保留字：数据结构面向目标的自适应a b st ra ct多物理场问题的成本和内存有效的有限元模拟仍然是一个具有挑战性的任务。由对偶加权残值法衍生出的面向目标的时空自适应方法是一项引人注目的关键技术。他们提供生成时空网格的最佳序列，以最大限度地降低成本。这项工作有助于有效地嵌入空时自适应方法的应用程序自由的方式数值筛选工具。本文结合DTM++项目提供的软件dwr-diffusion介绍了一种新的灵活的数据结构及其算法。©2019作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本1.0.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2018_262法律代码许可证https://github.com/dtm-project/dwr-diffusion/blob/master/LICENSE使用git的代码版本控制系统软件代码语言、工具和服务使用C++.17;cmake、gcc、mpi，可选：paraview、doxygen编译要求，操作环境依赖处理.II与hdf 5; Linux（Fedora，CentOS 7，RHEL 7等），MacOS，Windows WSL如果可用，链接到开发人员文档/手册https://github.com/dtm-project/dwr-diffusion问题支持电子邮件dtmproject@uwe.koecher.cc1. 动机和意义1.1. 介绍非均质变形多孔介质中流动问题的精确、可靠和高效的数值逼近这些现象包括多相扩散、对流主导的输运、化学反应以及开发高性价比的数值模拟筛选工具对支持环境、民用、能源、生物医学等众多工程领域的研究相当多的代表是基于不可压缩流的*通讯作者。电子邮件地址：koecher@hsu-hamburg.de（美国）。Köcher），bruchhaeuser@hsu-hamburg.de（M.P. Bruchhäuser），bause@hsu-hamburg.de（M.Bause）。https://doi.org/10.1016/j.softx.2019.100239耦合渗流与孔隙弹性波传播的Biot-Allard方程。Space–time [1，2]和参考文献。有限元法（FEM）可以网格化复杂的几何形状和边界，同时以自然的方式提供高阶近似。特别是间断Galerkin有限元法可以适应模型、材料参数和边界条件的间断性，比有限差分法和有限体积法更灵活。对于使用时空有限元方法的具有挑战性的问题，需要自适应和智能的网格选择来平衡精度和计算成本。面向目标的自适应性允许通过成本、能量或误差泛函而不是经典范数中的解变量来控制用户选择的数量中的误差。Becker和Rannacher介绍了双加权残差（DWR）方法，2352-7110/©2019作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx2联合 Köcher，M. P. Bruchhäuser和M. Bause / SoftwareX 10（2019）100239联系我们0∈∈0JJ−Jℓ=-= ≤∞n=τ，hn=Hn0000τ，h00j，0Ω0（）Q=×∈：Q→目标导向的适应性[3，4]。它被进一步应用和研究; cf.例如[5以及线性形式F：W→R，DWR方法的核心思想是将求解偏微分方程的问题嵌入到一个已知的问题中F（）：=联系我们0Ω·fdx dt，f∈L2（0，T;L2（英语：L2）））。最优控制的框架。目标是使成本、能量或感兴趣的误差函数最小化; Fig.1.该方法有助于找到例如最佳使整体计算成本最小化。它不简而言之，弱变分对偶问题是由下式给出的相应拉格朗日泛函的稳定条件引起的：依赖于初始时空网格的通常未知的假设。DWR方法具有挑战性的数学和计算的缺点相比，标准的后验误差J（u+）− J（u）− A（）（z）！F（n）−A（u）（n）{0}，（3）基于自适应方法。任意目标量的DWR误差估计在严格理论中不能保证下界的存在，这可能导致计算困难。此外，DWR方法需要解决两个问题，对于αV和z W，通过使用约束优化，参见[10，Sec. 6.1]、[9]。一般非线性目标泛函J（u）的对偶问题读作：找到z ∈ W，W=L2（0，T;H1（H）），满足z（T）=zT，zT∈H1（H），使得原始和辅助的更高的近似质量的对偶问题，导致更高的计算成本。这项工作提供了可扩展和有效的数据结构，A（z）（）=0J（u+）−J（u）−π0n（T）·z（T）dx，（4）一个通用的方式来结构和处理大量的计算数据，由于DWR的方法。DTM++.Project的相应软件dwr-diffusion实现了第1.2节给出的精确科学问题的独立求解器。第二部分介绍了软件体系结构和关键技术。在第3节中，详细介绍了该软件的一个示例。第4节解释了软件的影响，第5节通过结论完成了工作1.2. 软件解决的精确科学问题DTM ++.Project/dwr-diffusion frontend simulation tool forthe finite element analysis librarydeal.II[11]生成了非均质多孔介质中非稳态输运原型模型的参考实现。目标量是成本、误差或能量泛函。它取决于从原始解u导出的用户选择的感兴趣的量，即例如局部污染浓度或漂移。我们考虑目标量J（u）的目标导向近似，对任意的n∈V，V=L2（0，T;H1（n）），且伴随双线性形式A∈：W×V→R，A（z）（）：=·ρtz−·（ρz）dx dt，它是从A（z）的时间部分积分导出的，（2）并将其代入（3）的第一个等式。为了估计（u）（uτ，h）的函数误差并导出时空网格自适应性的可计算先验误差估计η原始问题和对偶问题通过适当的时空变分方法离散化，即，在这里，原始问题在时间上是分段常数不连续Galerkin方法，对偶问题在时间上是分段线性连续Petrov-Galerkin方法，与空间上的连续有限元方法相结合; cf. [9、10]。 时空圆柱CNOI被划分为不重叠Q = ×I，I<$= [t0，T]的分割为t0=：t <· ·<· ：=T，I ：=（t，t），t ：=t，n=1，. . .，N，for求u∈V 从扩散方程N的n mn mℓn−1第三圈。在每个Qn上，考虑a不一定符合分割T关于 分成非重叠元素K，表示为ρtu−·（u）=f，（1）h，n h n在有界），尺寸d=2（开放和如几何四边形对于d2或六面体对于d3。完全离散的原始和对偶解表示为：uτ，h，3和I（t0，T），0t0T，<并配备适当的初始和边界条件ℓτ，h（x，t）和z（x，t）。.循环次数=1，. . . 直到目标是u（x，t0）=u0（x）inn×{t0}，算法（dwr-扩散）到达; cf. 图1.一、解决原始问题：找到系数向量un，n =（u n，n）j，j = 1，. . . ，N primal，n，n，其中总数u（x，t）=g（x，t）在ΓD×I上，0j，0DOF原始的，n，n（x）=h（x，t）在ΓN×I上.rD<$rN=r N表示边界的划分。准确地说，空间N中的自由度（DoF）从DOF在板Q上，ΓD=表示Dirichlet部分，诺依曼部分（Mn+τnAn）un，n=τn（fn+hn）+Mn（Ih un（x，tn−1））j，边界为了简单起见，在下面的讨论中使用了关于ΓD=π的齐次Dirichlet边界值g=对于n1，. . . ，N，通过在时间上向前行进穿过石板。上原问题的质量和刚度矩阵软件实现包括非均匀边界这两种类型的值。系数函数ρ >0和ρ>0板n表示为Mn手侧组件fn和n和An，分别为右在（1）中，分别是介质的质量密度和渗透率，函数fR表示作用在内部区域上的体积力，例如源项或汇项。弱变分原问题读作：五，2 1 10h0 对 应 于 体 积 强 迫 和 非 齐 次Neumann边界项，并且向量（Ihu（x，tn-1））j是初始有限元上的初始值函数u0从（1）（n1）或前（n-1）个板（n>1）的当前板Q的空间。每个系统都经过修改，V=L（0，T;H0（H）），满足u（0）=u0，u0∈H0（H），使得A（u）（n）=F（n），n∈W=L2（0，T;H1（n）），（2）n强有力地应用来自（1）解决对偶问题：找到系数向量zn，n=（zn，n）j，j= 1，. . . ，N dual，n，n，总自由度具有双线性形式A：V×W→R，空间NDOF对偶，n，n对铸坯刘，从不Ω联合 Köcher，M. P. Bruchhäuser和M. Bause / SoftwareX 10（2019）1002393nn00nH τ，hTA（u）（）：=dxdt，0·ρDoFQn（2Md+τnAd）zn，τn =τnJn，τ n+2Md（Id zτ n（x，tn））Ω4联合 Köcher，M. P. Bruchhäuser和M. Bause / SoftwareX 10（2019）100239=nQ+=τ，h˜nτ，hh，n≤ ≤≤Fig. 1. 使用DWR技术进行误差估计η的dwr-扩散软件的实现摘要 （左）。一个模范的前进和后退时空片上的时间步进步长n （中图） 适应战略：空间三角测量更新 →T+1 时间三角测量前n精化I →I+1 并复制Q+1→Q+1 用于准备100 +1环（右）。h，nh，nn nnn+1对于n N，. . . ，1，通过在时间上穿越石板。板n上对偶问题的质量和刚度矩阵分别表示为Md和Ad1.4. 软件的基本设置和使用用户需要准备一个交易。IIv9.0 toolchain to com-n n侧装配向量Jn，堆并运行软件，这可以在任何主要平台上完成向量（Id z）0对应于用户选择的目标，并且（x，t n））j是初始值的插值通过使用candi;打开终端并键入Hτ，h函数zT（nN）或下一个（n1）板（nN）在当前板（n）的对偶有限元空间上的完全离散对偶解。<相同类型的边界着色（狄利克雷或诺依曼型）用于对偶问题，但具有齐次边界值函数，即使在非齐次原始边界条件的情况下。每个系统修改，使齐次Dirichlet边界条件git clonehttps://github.com/dealii/candicdcandi&&./ candi.sh并按照说明进行操作。然后下载、编译并运行gitclone https://github.org/dtm-project/dwr-diffusion cddwr-diffusion&&cmake.&&释放./ dwr-diffusion./输入/KoecherBruchhaeuser2d.prm对对偶解z∈D有很强的应用. 数字- n，n在第3节中给出的实验设置。cally近似的后验详细信息[9]和源代码文档。空间-时间网格细化更新如下实现：在时间上标记用于细化的空间-时间板n，其中η n，η n0≤θτ≤ 1的最大值，然后，在每个平板Q上标记一个网格1.5. 相关工作该代码的起源与用于拉普拉斯方程的deal.IIv8.4 [11]的第14步 教 程 代 码 相 关 并 从 中 派 生 。dwr-diffusion 软 件 使 用 其 他DTM++.Project（bitbucket.org/dtmproject）解算器的多种技术，胞腔K∈T<$对于空间中的d维各向同性精化，标量值，向量值和多物理问题作为介绍，其中，|η˜n,ℓ|属于顶部馏分第1.3节介绍。θ或θτ，h[11]的代码库提供了一些实现h，1h，2，对于没有或被标记用于时间细化的板，固定 弹塑性 问题 使用 基于DWR的目标-的最大值，0θh，2θh，11，然后执行空间细化，最后执行时间细化。1.3. 软件对科学发现的贡献本文介绍的dwr扩散软件是DTM++项目的一部分。DTM++.Project为deal.II有限元库收集前端解算器。该集合提供了高阶时间离散化、复杂求解器和预处理技术的实现，以模块化的方式支持用于非定常PDE问题的分布式存储器并行的参数输入和数据输出的类几个熟悉的DTM++的前端解算器。项目套件，例如。用于声波、弹性波和耦合波传播的XWAVE，用于质量守恒传输的MEAT和Biot耦合变形与传输，已经对科学发现的过程做出了贡献;[12]和参考文献。[14 ]中使用了dwr扩散的固定前趋，[ 9 ]中使用了后继，用于稳定的非固定对流主导的传输问题。贡献的软件共享所有DTM++的模块化和灵活性。项目求解器，它提供了一个免费可用的开源框架，用于有效地解决具有目标导向适应性的问题，并将支持进一步的科学发现与物理相关的挑战性问题的数值模拟。定向网格自适应。它还提供了非定常对流扩散反应方程的实现没有自 适应 。另一个基于deal.II 的模 块化 开源软 件框架 是DOpElib（dopelib.net），它为（当前）静态问题提供了DWR方法，特别强调解决优化问题。虽然DOpElib提供了一个具有最少用户输入的框架来解决PDE或具有PDE约束的优化问题，但人们不能轻易地将数据结构从该框架中提取到独立的软件中。所有列出的相关实现都可以从我们的工作中获益，因为这些实现大多是特定的实现，部分缺乏文档，更重要的是，没有以无应用程序的方式记录它们的数据结构和算法2. 软件描述DTM++.Project的软件dwr-diffusion以面向对象的方式实现了第1.2该软件是用C++.17语言编写的，并利用了C++.11以来的新语言功能，如引用计数动态内存分配，基于范围的循环，自动说明符，强类型枚举类等。联合 Köcher，M. P. Bruchhäuser和M. Bause / SoftwareX 10（2019）1002395J+τ，hQτ，Q= × Ωnh，1nτ，h图二. dwr-扩散数据存储管理的图示概述。列表条目对应于slab，并存储共享指针引用的小数组。大数据Vector独立存储，可以串行化。2.1. 软件构架dwr-diffusion软件是deal.IIv9.0库的前端求解器[11]。这是7个例子。这些示例是通过从命令行调用相应的基于文本的参数输入文件来驱动的。运的例子来自给定的解析解，以确保实现的正确性。第 1.2 节 中 的 算 法 的 实 现 在 单 个 中 心 求 解 器 类 模 板Diffusion_DWR dim>中完成。这个基类是从DTM：：Problem继承的，用于实现无量纲方法。这种方法可以扩展到支持不同的时间积分方法;参见。[12 ]第10段。求解器类的结构模仿算法给出的顺序，图1.一、给出了原始和对偶问题的时间离散化略有不同的统一实现。准确地说，时间基函数的求积规则是图二、通过一种轻量级的数据结构访问大量的空时自由度矢量数据。数据存储man结构只需要存储引用计数指针。列表数据结构允许有效地向前和向后搜索以访问相邻空时片上的解决方案这对于给定的向前或向后时间推进实现是需要的未来的发展是一个智能存储的man-serialises类的时空自由度矢量的权重大的数据，只要他们不需要实际的计算。我们在这里不提供这样的类，因为这将是一个特定于应用程序的决定。我们的数据存储方法是，甚至然后（内存）高效的孔时空原始解决方案u是需要的建设的后一种情况可能发生在非线性情况下; cf. [5] Rem. 3.1]。另一个重要的关键技术是新的DTM：：Grid_DWR类的列表方法，类似于数据存储管理。 此类允许处理（d1）维时空计算板。它允许廉价的本地时间优化。对于每个板，需要有效地存储和访问空间和时间三角测量、有限元对象、约束对象和DoF处理程序;参见图10。图3和1.出于效率原因，不存储诸如系统矩阵组件的其他数据3. 说明性实例软件实现目标函数J：V→R，J（）：=·（u−u ）dx dt，（5）可以在参数输入的高级别选择生成的参数。所有的汇编器，例如用于矩阵的，右边的vec-乌乌河中国cIc布拉奇τ，htors甚至对于所有的v∈V，u∈ v，表示第二个问题用于线程并行和进一步MPI+X并行模拟的汇编器被外包到单独的类中，以允许它们在其他框架中重用，只需最少的编码工作。因此，DTM++.Project套件的有用和有效的辅助类是独立提供的，因此输入参数处理和数据输出处理不必包含在求解器类中。PDE求解器框架的通用DTM++.Project软件框架概述可参见[12]和代码的源代码doxygen文档DWR循环，u表示给定的解析解，其目的是达到u−u<对于公差标准tol，在控制体积c上I c;cf. 图 四、备注。提供了一个类似范数的目标泛函，允许一个简单的比较面向目标的误差估计与其他误差估计。为了证明这些功能，一个解析解u，它模拟了一个具有随时间变化的高度的旋转圆锥u（x，t）：=u1·u2，x=[x1，x2]∈R2且t∈R，2.2. 软件关键技术u1（x，t）：=（1+a·（（x1−m1（t））2+（x2−m2（t））2））−1，u2（t）：= ν1（t）·s·arctan（ν2（t）），（七）该工作的关键技术之一是有效地处理时空自由度数据的近似于二维L形域（图4），I =（0，1. KB2D测试用例。辅助数据图3.第三章。 处理计算板的 元 素的 图 示 概 述 Q= T 通过dw r-d if fu s ion 软 件 的DTM：：Grid_DWR类实现的I类。6联合 Köcher，M. P. Bruchhäuser和M. Bause / SoftwareX 10（2019）100239n=×}=\τ，=·−|η˜τ，h∈≈τ，2τ，42412h，h，τ，τ，（n个1）其中，u− u18 = 5。744· 10−4。见图4。实施的局部（d+1）维<时空控制体积Qc（t）= Ic（t）× Ic <$Q =Ic × I的图示，其中Ic =（t m，c，t n，c）并且对于等式t 0 ≤ t m，c t n，c≤ T。（五）第三节。图六、第3 节时 空片Q <$ ，<$ ∈ { 17 ，18 }的τn分布。定义为m1（t）：=1+1cos（2πt）和m2（t）：=1+1sin（2πt），并且，ν1（t）：=−1，ν2（t）：=5π·（4t−1），其中t∈[0，0.5）和ν（t）：=1，ν（t）：=5π·（4（t−0. 5）−1），对于t∈[0.5，1），t=t−k，表1选择的收敛性和误差估计结果的环路滤波器，具有N个slab和最大。 网格单元数 |T|对于第3节，在平板n = 100上。k∈ N0，标量a，s ∈ R，a> 0. 一个依赖时间的空间吴恩达|T|u−u时间控制量Qc，KB2d∈c（t）Ic。c（t）由下式给出c（t）={x∈R2|x−[0. 五，零。[5]T=x<$+m（t;r1），x<$∈m，m（t;r1）=r1·[cos（t），sin（t）]T}，1 5 3 6.070e−02 1.994e−02 0.332 6 12 2.634e− 02 2.451e−02 0.931·72274417··3·04e−042.030e−040.2818 295 603 5.744e−04 1.723e−04 0.30式中，I_c=（p_（1，2））×（p_（1，2））R_（1，2），I_c=（0. （见第25页，第1页）。1 1 2 2224. 影响该设置使用了{（0，1）}\{（0. 5，1）}，ΓN={x ∈H|X1=0，ΓDhIN. 导出了（1）中的函数u0，f，g和h根据解析解（7），将系数函数设置为ρ= 0。8，则x= 1。2，标量设为a= 50，s=−10。3333. 构造了φ = 1的的N =5个平板，每个平板具有3个空间单元。目标公差是tol10−2uu1h Qc. uτ、h和等高线的解曲线zτ，h的线如图所示。五、选定的收敛和误差估计结果以及质量测量Ieff：=软件dwr-diffusion软件为目标导向网格自适应性的实现提供了新的有效见解，支持基于有限元方法的经济高效的数值筛选工具的开发，以解决科学和工程中的基本问题。该软件的衍生物已经用于稳定的固定和非固定对流主导的问题，ℓτ，h/（J（u）−J（u）））|在那里，我 =1是最优的，给出 通过网格自适应性[9，14]以及耦合变形和扩散驱动的运输[12，13]。的目标用户群表1. 图6示出了Q的时间子区间长度在I上的分布。当τ=18时，达到了判据（6）数值结果（Fig. 5、突出地方特色--以自适应方式优化求解器的时间细化特性，而更精细的时空单元位于控制体积Qc，KB2d（x，t），t（0. （见第25页，第 1页）。请注意，由于控制体积不活动，因此对于该抛物线问题的t>1不进行任何细化备注。KB2d测试用例（表1）的质量测量Ieff似乎不太令人信服。这是由于我们已经针对具有与时间无关的控制体积的更简单的情况KB2、KB2a、KB2b调整了（静态）参数θτ、θh，1、θh，2以达到Ieff 1。对于KB2d情况，控制体积随时间变化，不容易找到最佳静态调谐参数。 我们应该在这里优化θτ（n），θh，1（n，n），θh，2（n，n），并且可以对每个n和n取不同的值以达到Ieff= 1。0，这超出了本工作的范围。研究人员和工程师使用和开发基于有限元的数值模拟筛选工具来解决几个问题是广泛的和全球分布的;一个例子是工作[15]的扩展，其中他们使用高阶变分时间离散化来进行对流主导的传输，具有非目标导向的自适应时间步长控制，或者是对[7]中提出的方法的有效非平稳扩展。该软件目前未用于商业环境，但许可证允许集成到商业工具中。5. 结论DTM++.Project的原始软件出版物dwr-diffusion以高度模块化的方式为实现面向目标的网格自适应性提供了高效和可扩展的数据结构和算法。数据结构和算法的关键思想可以被重用，甚至可以独立地重用 在任何自适应有限元图五. 解决方案配置文件u18和双解z18轮廓线的第3部分。联合 Köcher，M. P. Bruchhäuser和M. Bause / SoftwareX 10（2019）1002397代码. 该软件的性能和适用性显示了一个说明性的例子和其他几个人的代码。这项工作产生了一个重大的突破，作为一个免费提供的开源实现与广泛的源代码文档的双重加权残差方法的insta-tionary问题，在一个应用程序免费的方式，使它可以很容易地被其他框架采用。正在进行的工作是分布式内存并行化，隐藏内存操作的数据存储的串行化，并找到最佳的调整参数的最小化问题。附录A. 补充数据与本文相关的补充材料可以在https://doi.org/10.1016/j.softx.2019.100239上找到。引用[1]Hughes TJR，Hulbert GM.弹性动力学的时空有限元法：公式和误差估计。计算方法应用机 械 工 程 1988;66（3）：339-63.[2]放大图片作者：Schwab C.抛物型发展问题的时空自适应小波方法。MathComp2009;78（267）：1293-318.[3]贝克河，兰纳赫河有限元方法中的加权后验误差控制。In：Bock HG，et al.，编辑器。第二届EC会议的程序编号。数学和Adv. Appl. ENUMATH 1997.新加坡：世界科学; 1998年。[4]贝 克 河 ， 兰 纳 赫 河 有 限 元 法 中 后 验 误 差 估 计 的 最 优 控 制 方 法 。 ActaNumer2001;10：1-102.[5]Besier M，Rannacher R. 计算非定常不可压缩流的有限元Galerkin方法中面向目标的时空自适应性。Int J Numer Methods Fluids2012;70（9）：1139-66.[6]里希特·迈锡尼不可压Navier-Stokes方程分数步θ离散的后验误差估计。计算方法应用机械工程2015;288：45-59.[7]Endtmayer B，Wick T.椭圆型问题多目标误差估计的单位分解双加权残值法。ComputMethods Appl Math2017;17（4）：1-25.[8]作者：Wick T.非线性流固耦合的自适应时步控制。J Comput Phys2018;366：448-77.[9]放大图片作者：Bruchhäuser MP，Schwegler K，Bause M.非定常对流占优方程基于双加权残差的误差控制：势还是压载？。2018，p. 1[10]作者：J.微分方程的自适应有限元方法。Basel：Birkhäuser;2003.[11] Alzetta G，Arndt D，Bangerth W，Boddu V，Brands B，Davydov D，Gassmoeller R ， Heister T ， Heltai L ， Kormann K ， Kronbichler M ，Maier M，Pelteret J-P，Turcksin B，Wells D. 协议。II库，9.0版。 JNumerMath 2018;26（4）：173-83.[12]克歇尔大学弹性波方程和扩散方程的变分时空方法（Ph.D. thesis），Mech.Engrg. Helmut-Schmidt-University Hamburg; 2015，p. 1-188[13]Bause M，Radu FA，Köcher U. Biot多孔弹性系统迭代耦合的时空有限元逼近。计算方法应用机械工程2017;320：745-68。[14]放大图片作者：Bruchhäuser MP，Schwegler K，Bause M.稳定化有限元法目标导向误差控制之数值研究。In：Apel T，et al.， 编辑器。高级有限元法与计算机应用讲义。Sci. 和工程师，Berlin：Springer; 2018，p. 1[15]Ahmed N，John V.高阶变分时间离散应用于对流扩散方程的自适应时间步长控制。计算方法应用机械工程2015;285：83-101.