面向高阶精度CFD的JFNK算法及其并行计算*

程 彬，李大力，徐传福，刘 巍，王光学，邓小刚

1.国防科学技术大学 计算机学院，长沙 410073

2.中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000

3.国防科学技术大学，长沙 410073

面向高阶精度CFD的JFNK算法及其并行计算*

程 彬1+，李大力1，徐传福1，刘 巍1，王光学2，邓小刚3

1.国防科学技术大学 计算机学院，长沙 410073

2.中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000

3.国防科学技术大学，长沙 410073

CHENG Bin,LI Dali,XU Chuanfu,et al.Research on Jacobian-Free Newton-Krylov method for high-order CFD applications and its parallel computing.Journal of Frontiers of Computer Science and Technology, 2017,11(1)：61-69.

目前计算效率低是限制计算流体力学（computational fluid dynamics，CFD）高阶精度格式方法的重要因素之一。由于高阶精度格式计算模板相对复杂，很难精确计算其Jacobian矩阵，从而影响传统LU-SGS（lowerupper symmetric Gauss-Seidel）等算法的收敛效率。JFNK（Jacobian-free Newton-Krylov）算法是Krylov子空间方法与非精确牛顿方法的结合，拥有较好的迭代收敛效率，采用无矩阵思想，只计算Jacobian矩阵与矢量的乘积，从而有效避免Jacobian矩阵的计算和存储。在真实高精度结构网格CFD应用程序中，设计并实现了JFNK时间求解算法。在有粘低速圆柱绕流的算例测试中，和传统LU-SGS算法相比，JFNK算法拥有更好的计算稳定性，同时可使迭代收敛效率提高2倍以上。以天河2号超级计算机为并行计算平台，对JFNK算法和传统的LU-SGS算法的并行强可扩展性进行了测试，二者均表现出良好的并行效率。

高阶精度格式；JFNK算法；计算效率；可扩展性

1 引言

计算流体力学（computationalfluid dynamics，CFD）是一门采用数值方法求解流动控制方程，从而揭示流动本质和规律的学科[1]。随着CFD的发展以及高阶精度格式的应用，可模拟的流动问题越来越复杂，时空尺度越来越精细，迫切需要提高模拟效率以加快复杂问题的模拟过程，这就对时间求解算法和并行计算提出了更高的要求。

高精度格式具有计算精度高，流场分辨率高的优点，能够很好地适用于激波、湍流、边界层等复杂流体问题的模拟。Harten等人提出了本质无震荡格式[2]（essentially non-oscillatory scheme，ENO）；Liu等人基于ENO格式提出了加权本质无震荡格式[3]weighted essentially non-oscillatory scheme，WENO）；邓小刚等人提出了加权紧致非线性格式[4-5]，包括WCNS（weighted compact nonlinear scheme）、HWCNS（hybrid cell-edge and cell-node weighted compact nonlinear scheme）以及混合型耗散紧致格式[6]hybrid cell-edge and cell-node dissipative compact scheme，HDCS）；Reed和Hill提出了间断有限元方法[7]discontinuous Galerkin method，DG）。但是相对于二阶精度格式，高阶精度格式单时间步计算量较大，计算稳定性较差，导致高阶精度格式的模拟效率较低，甚至不容易收敛。解决这一问题则需要从时间求解算法和并行计算两个方面来综合考虑。

在算法方面，采用高效的迭代求解算法能够有效减少迭代步数，从而提高CFD模拟的效率，经典的时间求解算法有点松弛SGS（point relaxation symmetric Gauss-Seidel）、线松弛SGS（line relaxation symmetric Gauss-Seidel）、LU-SGS（lower-upper symmetric Gauss-Seidel）、LU-ADI（lower-upper alternating direction implicit）和Jacobi迭代法。经典迭代法计算稳定性高，应用广泛，然而随着计算规模持续膨胀，经典迭代法计算效率受到限制，无法满足工程实际问题需求。Krylov子空间方法作为一种高效的线性系统求解方法在近年来得到了广泛的研究和应用。Krylov子空间方法是对基于Krylov子空间的一类算法的总称，包含很多种不同的解法，常用的有CG（conjugate gra-dient）、GMRES（generatized minimal residual）、JFNK（Jacobian-free Newton-Krylov）等。其中，JFNK算法是采用非线性近似计算Jacobian矩阵和矢量乘积的Newton-Krylov子空间方法。研究表明预条件能够有效提高JFNK算法的收敛速度[8-9]。预条件技术将方程组变换为另一个更容易迭代求解的同解方程组。其实质就是在保持方程组解不变的前提下，缩小矩阵特征值范围，以改善矩阵的性质，从而达到加速收敛的目的。McHugh等人[10-12]研究了低阶空间离散格式CFD应用中各种预条件子在Krylov子空间方法中的效率；Knoll等人[13]研究了多重网格法作为预条件子，将JFNK算法应用于求解不可压Navier-Stokes方程。上述方法都得到了较好的加速收敛效果。不过关于Krylov子空间方法及其预条件在高阶精度格式中的收敛性及稳定性研究还不够深入全面。

CFD并行计算的基本思想是采用区域分解的数据并行方法，把流场区域划分为不同的子区域，将流场信息、几何信息载入相应区域的CPU内存中，每个子区域分配到节点上计算，在子区域内进行串行迭代，一步计算完成后子区域之间交换边界，然后进行下一步计算，直到收敛。并行算法的并行扩展性对CFD模拟效率同样有着重要的影响。

本文基于自主开发的真实复杂结构网格高阶精度CFD程序HOSTA（high-order simulator for aerodynamics）研究了JFNK求解方法及其并行计算技术。HOSTA采用有限差分法离散控制方程，空间离散采用WCNS[14]高阶精度格式，在多块结构网格上求解Navier-Stokes方程，已应用于航空航天领域等复杂外形复杂流动问题的机理研究和工程应用中。HOSTA中实现了一系列的定常、非定常时间求解方法，包括三步Runge-Kutta、LU-SGS、标量点松弛、矩阵点松弛等。本文在原有时间求解方法上实现了预条件JFNK求解方法，采用圆柱绕流算例深入研究分析了JFNK算法在高精度格式WCNS中的迭代收敛效率和算法稳定性，同时在天河2号超级计算机上实现了JFNK求解算法的并行计算，分析了HOSTA程序的计算效率和并行可扩展性。

2 高阶精度格式常用时间求解算法

计算坐标系下，不考虑源项，无量纲Navier-Stokes方程表示为[15]：

式中，为无量纲形式的守恒变量矢量；是l(l=ξ,η,ζ)方向的无粘通量矢量；是l(l=ξ,η,ζ)方向的粘性通量矢量。由式（1）可得：

式中，J为坐标变换Jacobian行列式。

时间求解的收敛速度由左手项决定，一般采用Jacobi、LU-ADI、Gauss-Seidel等线性方程组方法或者Krylov子空间方法求解。其中LU-SGS隐式格式是Navier-Stokes方程定常问题中隐式时间推进方法应用最广泛的[16]。LU-SGS方法构造简单，在二阶格式中具有稳定性好和收敛效率高的特点；但在高阶格式中，由于右端项采用高阶精度格式离散，左端项和右端项离散精度不匹配，影响其收敛性和收敛效率[17]。

3 JFNK算法及预条件

JFNK算法是求解大型非线性方程组的算法[10,18]，是非精确Newton法与Krylov子空间法的结合。JFNK算法是在GMRES算法的基础上，采用无矩阵思想，利用有限差商计算来近似Jacobian矩阵和矢量的乘积计算，因此不需要计算和存储Jacobian矩阵。在时间步迭代之内，JFNK算法至少包含两层迭代，外层是非精确Newton迭代，内层是Krylov子空间线性迭代。本文着眼于JFNK算法的预条件等参数对于计算效率的影响，同时对比分析JFNK算法和LU-SGS算法的并行扩展性。

3.1 JFNK算法

JFNK算法的计算步骤如下：

式中，m代表此算法最大内层迭代数；ε是最小范数收敛判据；r0=b-Ax0和wj=Avj采用无矩阵思想，利用有限差商计算即为JFNK算法。

3.2 Jacobian矩阵和矢量乘积的计算

由于Krylov子空间方法只用到Jacobian矩阵与向量的乘积，从而可以不计算和存储Jacobian矩阵。Krylov空间记为：

令v=Ai-1r(k)，则可采用有限差商近似计算Jacobian矩阵与向量的乘积：

式中，ε对于方程求解的精度和稳定性影响很大，一种计算方法[19]是：

3.3 JFNK算法的预条件

对于复杂数值计算，JFNK算法迭代中需要占用大量计算资源，为了加快迭代收敛速度，考虑到特征值分布越集中收敛越快的特点，在求解过程中，可以采用预条件来改善线性系统。

例如，引入线性左预条件，写成典型线性方程组形式，然后左乘M-1得：

式中，M是矩阵A的近似矩阵，M-1位于线性层次，迭代过程中不变，为M(Q(k))-1。预条件后的方程残差为：

预条件后的Krylov子空间为：

Km(M-1A,r(k))≡span{r(k),M-1Ar(k),…,(M-1A)m-1r(k)}采用JFNK算法时，令v=(M-1A)i-1r(k)，则：

本文采用了标量和矩阵的LU-SGS作为JFNK算法的预条件，二者在求解预条件方程时对Jacobian近似矩阵的分裂上存在差异。矩阵预条件采用的是按特征值分裂，而标量预条件采用的是按谱半径分裂。关于标量和矩阵方法的介绍见文献[15]。

4 JFNK算法及预条件的并行实现

并行计算方面，HOSTA程序采用了区域分解的思想结合MPI（message passing interface）并行编程实现多进程的并行计算，每次迭代计算过程中仅对边界信息进行通信交互。

当引入JFNK算法后，需要采用基向量构造Krylov子空间，因此在JFNK算法之前要将原来的多维数据结构映射到一维向量中。采用各块按序分段映射的方法，一个网格块对应向量中的一段，保证向量的分段与区域分解的分块对应，同时又方便块内的内积运算。

如3.1节的算法流程中所述，JFNK算法中存在向量内积运算，因此需要引入进程归约通信。首先各个进程对所拥有的向量部分做内积运算；然后通过MPI归约操作，将各个进程上的局部内积累加到根进程；最后再通过MPI广播操作将最终的内积结果广播给所有进程。

另外如图1和图2所示，预条件的引入还将涉及到预条件方程右端项的计算（包括通信）以及预条件方程的求解等过程。

5 算例与测试平台介绍

本文所采用的测试算例是定常低速圆柱绕流问题，分析对比预条件JFNK算法和LU-SGS算法的收敛效率及其影响因素。采用的是单块网格，网格规模为96×104×5（径向×周向×展向），网格截面图如图3所示。圆柱的无量纲直径和高度等于1，第一层无量纲网格高度等于0.05。计算状态为马赫数Ma=0.1，雷诺数Re=5，T=287 K。使用3种隐式时间推进算法，包括标量LU-SGS算法以及标量预条件和矩阵预条件的JFNK算法。图4为收敛流场无量纲速度U的分布图。本文还对算法的强并行可扩展性进行了测试分析，算例网格规模为512×512×5（径向×周向×展向），利用区域分解技术将网格均匀分割。

Fig.1 Flow chart of HOSTAprogram图1 HOSTA程序流程图

Fig.2 Flow chart of HOSTAprogram using JFNK algorithm图2 使用JFNK算法的HOSTA程序流程图

Fig.3 Computational grid of flow around a circular cylinder图3 圆柱绕流算例计算网格

Fig.4 Dimensionless velocity contour distribution of flow around a circular cylinder图4 圆柱绕流算例无量纲速度分布

本文的测试平台是天河2号超级计算机系统。天河2系统由16 000个异构体系结构的计算节点组成，总处理器核数为312万个，双精度浮点运算的总峰值性能为54 PFlop/s，2013年6月至今一直位居全球超级计算机500强排名榜首。每个天河2号计算节点由2块Intel Xeon E5-2692处理器和3块Intel Xeon Phi 31S1P（MIC）协处理器组成。本文测试仅在CPU同构平台进行，每个计算节点有24个CPU核，共享64GB内存空间，双精度峰值浮点性能高达422 GFlop/s。

6 收敛效率测试分析

图5和图6分别给出了平均残差相对于时间迭代步和墙钟时间的平均残差收敛曲线，结果表明不管是在迭代步数还是墙钟时间上JFNK算法都要优于LU-SGS算法。JFNK算法残差达到收敛的时间步数仅为LU-SGS算法的1/4，墙钟时间也缩短了超过一半，其中矩阵预条件JFNK算法迭代步和墙钟时间略多于标量预条件JFNK算法。

Fig.5 Average residual curves vs.time step number of LU-SGS and JFNK with nonlinear scalar and matrix type preconditioner图5 标量、矩阵预条件JFNK和LU-SGS关于时间步的平均残差曲线

Fig.6 Average residual curves vs.wall time of LU-SGS and JFNK with nonlinear scalar and matrix type preconditioner图6 标量、矩阵预条件JFNK和LU-SGS关于墙钟时间的平均残差曲线

图7和图8给出了CFL数对标量预条件和矩阵预条件JFNK算法收敛效率的影响，总体上CFL数越大算法的收敛性越好。当CFL数较小时，迭代收敛速度随CFL数变化较快，当CFL数较大时，迭代收敛速度随CFL数变化较慢，最终，收敛速度基本保持不变。对标量预条件和矩阵预条件JFNK算法测试都体现了类似的规律。同时还发现JFNK算法可以取到比LU-SGS算法更大的CFL数，即JFNK算法体现出更优的计算稳定性。综上分析，JFNK算法拥有相对LU-SGS算法更好的收敛效率和计算稳定性。

Fig.7 Influence of CFL number on convergence rate of JFNK with nonlinear scalar type preconditioner图7 CFL数对标量预条件JFNK收敛速度的影响

Fig.8 Influence of CFL number on convergence rate of JFNK with nonlinear matrix type preconditioner图8 CFL数对矩阵预条件JFNK收敛速度的影响

7 并行扩展性分析

图9给出了单节点内不同进程数的程序强可扩展性结果，其中网格分为24块，在单节点上采用如横坐标所示的进程数对其进行测试。从图中可以看出，在进程数小于4时，并行效率可以保持在80%以上，进程数大于8时，受访存带宽限制，程序的并行效率下降明显。JFNK算法和LU-SGS算法具有类似的问题，因此在进行多节点测试时，每个节点的进程数最多取到4。

Fig.9 Strong program scalability on a single node图9 单个节点强可扩展性测试

图10给出了多个节点不同进程的程序强可扩展性结果，其中网格分为64块，在16个节点上对进程扩展性进行测试，最大进程数为64。从图中可以看出，随着进程数增加，程序的并行效率能保持在85%以上，JFNK算法和LU-SGS算法的并行效率非常接近，都具有较好的强可扩展性。

Fig.10 Strong program scalability on multiple nodes图10 多个节点强可扩展性测试

8 结束语

高阶精度格式CFD有限差分模型的计算是非常复杂的，由于难以得到精确的Jacobian矩阵，采用近似Jacobian矩阵的LU-SGS算法收敛效率不高，也无法与高阶空间离散格式匹配。JFNK算法避免了计算和存储Jacobian矩阵，不仅收敛效率更高，同时表现出了更好的稳定性，尤其适合于采用高阶空间离散的CFD应用。本文测试算例表明，JFNK算法可以取到比LU-SGS算法更大的CFL数，使用标量预条件、矩阵预条件的JFNK算法都要好于LU-SGS算法，墙钟时间节省超过一半。

在天河2号超级计算机上进行了并行测试，结果表明JFNK算法和LU-SGS算法都具有良好的强可扩展性，并行计算可以有效地提高HOSTA程序的计算效率，HOSTA程序在高性能应用方面具有优势。

后续工作除了进一步分析JFNK算法的其他影响因素，包括内层迭代数、不同类型的预处理条件、预处理迭代数、松弛因子等对高阶精度CFD应用收敛效率和稳定性的影响，同时将在天河2号超级计算机上对JFNK算法实现CPU+MIC异构协同并行计算和深度性能优化。

[1]Yan Chao.Methodology and applications of computational fluid fynamics[M].Beijing:Beijing University of Aeronautics andAeronautics Press,2006.

[2]Harten A,Engquist B,Osher S,et al.Uniformly high order accurate essentially non-oscillatory schemes[J].Journal of Computational Physics,1987,71(2):231-303.

[3]Liu Xudong,Osher S,Chan T.Weighted essentially nonoscillatory schemes[J].Journal of Computational Physics, 1994,115(1):200-212.

[4]Deng Xiaogang,Mao Meiliang.Weighted compact high-order nonlinear schemes for the Euler equations[C]//Proceedings of the 13th AIAA Computational Fluid Dynamics Conference,Snowmass Village,USA,1997.Reston,USA:AIAA, 1997.

[5]Deng Xiaogang,Mao Meiliang,Jiang Yi,et al.New highorder hybrid cell-edge and cell-node weighted compact nonlinear schemes[C]//Proceedings of the 20th AIAA Computational Fluid Dynamics Conference,Honolulu,USA,Jun 27-30,2011.Reston,USA:AIAA,2011.

[6]Deng Xiaogang,Jiang Yi,Mao Meiliang,et al.Developing hybrid cell-edge and cell-node dissipative compact scheme for complex geometry flows[J].Science China Technological Sciences,2013,56(10):2361-2369.

[7]Zhang Laiping,Liu Wei,He Lixin,et al.A class of hybrid dg/fv methods for conservation laws I:basic formulation and one-dimensional systems[J].Journal of Computational Physics,2012,231(4):1081-1103.

[8]Wu Jianping,Wang Zhenghua,Li Xiaomei.Efficient solving and parallel computing of sparse linear systems[M].Changsha:Hunan Science and Technology Press,2004.

[9]Cai Xiaochuan,Keyes D E.Nonlinearly preconditioned inexact Newton algorithms[J].SIAM Journal on Scientific Computing,2001,24(1):183-200.

[10]Johnson R W,McHugh P R,Knoll D A.High-order scheme implementation using Newton-Krylov solution methods[J]. Numerical Heat Transfer:Part B Fundamentals,1997,31 (3):295-312.

[11]McHugh P R,Knoll D A.Fully implicit finite volume solutions of the incompressible Navier-Stokes and energy equations using inexact Newton's method[J].International Journal for Numerical Methods in Fluids,1994,19(5):439-455. [12]McHugh Knoll P R,D A.Inexact Newton's method solution to the incompressible Navier-Stokes and energy equations using standard and matrix-free implementations[C]// Proceedings of the 11th Computational Fluid Dynamics Conference,Orlando,USA,1994.Reston,USA:AIAA,1994.

[13]Knoll D A,Rider W J.A multigrid preconditioned Newton-Krylov method[J].SIAM Journal on Scientific Computing, 1999,21(2):691-710.

[14]Che Yonggang,Zhang Lilun,Wang Yongxian,et al.OpenMP performance analysis of CFD application on Intel multicore and manycore architectures[J].Journal of Frontiers of Computer Science and Technology,2015,9(10):1153-1162.

[15]Liu Wei,Zhang Lilun,Wang Yongxian,et al.Foundations of computational serodynamics parallel programming[M]. Beijing:National Defense Industry Press,2013.

[16]Jameson A,Yoon S.Lower-upper implicit schemes with multiple grids for the Euler equations[J].AIAA Journal, 1987,25(7):929-935.

[17]Klopfer G H,Hung C M,van der Wijgaart R F,et al.A diagonalized diagonal dominant alternating direction implicit (D3ADI)scheme and subiteration correction[C]//Proceedings of the 29th AIAA Fluid Dynamics Conference,Albuquerque,USA,1998.Reston,USA:AIAA,1998.

[18]Brown P N,Saad Y.Hybrid Krylov methods for nonlinear systems of equations[J].SIAM Journal on Scientific and Statistical Computing,1990,11(3):450-481.

[19]Park H K,Nourgaliev R R,Martineau R C,et al.On physicsbased preconditioning of the Navier-Stokes equations[J]. Journal of Computational Physics,2009,228(24):9131-9146.

附中文参考文献：

[1]阎超.计算流体力学方法及应用[M].北京:北京航空航天大学出版社,2006.

[8]吴建平,王正华,李晓梅.稀疏线性方程组的高效求解与并行计算[M].长沙:湖南科学技术出版社,2004.

[14]车永刚,张理论,王勇献,等.Intel多核与集成众核上CFD程序的OpenMP性能分析[J].计算机科学与探索,2015,9 (10):1153-1162.

[15]刘巍,张理论,王勇献,等.计算空气动力学并行编程基础[M].北京:国防工业出版社,2013.

CHENG Bin was born in 1991.He received the M.S.degree in computer science from National University of Defense Technology.His research interest is applications of high performance computing.

程彬（1991—），男，湖南澧县人，国防科技大学计算机科学与技术专业硕士，主要研究领域为高性能计算应用。

LI Dali was born in 1987.He is a Ph.D.candidate in computer science at National University of Defense Technology. His research interest is applications of high performance computing.

李大力（1987—），男，湖北荆州人，国防科技大学计算机科学与技术专业博士研究生，主要研究领域为高性能计算应用。

XU Chuanfu was born in 1980.He received the Ph.D.degree in computer science from National University of Defense Technology.Now he is an assistant professor at College of Computer,National University of Defense Technology.His research interests include parallel algorithm and program optimization.

徐传福（1980—），男，安徽六安人，国防科技大学计算机科学与技术专业博士，现为国防科技大学计算机学院助理研究员，主要研究领域为并行算法，程序优化。

LIU Wei was born in 1980.He received the Ph.D.degree in aeronautical and astronautical science and technology from National University of Defense Technology.Now he is an assistant professor at College of Computer,National University of Defense Technology.His research interest is computational fluid dynamics.

刘巍（1980—），男，吉林辽源人，国防科技大学航空宇航科学与技术专业博士，现为国防科技大学计算机学院助理研究员，主要研究领域为计算流体力学。

WANG Guangxue was born in 1976.He is an associate professor at State Key Laboratory of Aerodynamics,China Aerodynamics Research and Development Center.His research interest is computational fluid dynamics.

王光学（1976—），男，重庆人，中国空气动力研究与发展中心空气动力学国家重点实验室副研究员，主要研究领域为计算流体力学。

DENG Xiaogang was born in 1960.He is a professor and Ph.D.supervisor at National University of Defense Technology.His research interest is computational fluid dynamics.

邓小刚（1960—），男，四川绵阳人，国防科技大学教授、博士生导师，主要研究领域为计算流体力学。

Research on Jacobian-Free Newton-Krylov Method for High-Order CFD Applications and Its Parallel Computing*

CHENG Bin1+,LI Dali1,XU Chuanfu1,LIU Wei1,WANG Guangxue2,DENG Xiaogang3
1.College of Computer,National University of Defense Technology,Changsha 410073,China
2.State Key Laboratory of Aerodynamics,China Aerodynamics Research and Development Center,Mianyang,Sichuan 621000,China
3.National University of Defense Technology,Changsha 410073,China
+Corresponding author:E-mail:13297499104@163.com

Computational efficiency is a significant factor of high-order accurate scheme for computational fluid dynamics(CFD)finite difference method.It is hard to get Jacobian matrix for lower-upper symmetric Gauss-Seidel(LUSGS)method because of its complicated computing stencil,which affects the computational efficiency.Jacobian-free Newton-Krylov(JFNK)method is a combination of inexact Newton method and Krylov subspace method.With the adoption of matrix free technique,JFNK method only calculates the approximate product of Jacobian matrix and vec-tor using the finite difference quotient of the nonlinear function,thus successfully avoids the calculation and store of approximate Jacobian matrix.Based on the domestic high-order simulator for aerodynamics program,this paper designs and realizes JFNK method.In the viscid and low speed test case of steady flow around a circular cylinder,compared with LU-SGS time stepping method,JFNK method is more stable and efficient.Some numerical experiments of JFNK method and LU-SGS method are performed on Tianhe-2 supercomputer system and good scalability is observed from the test results.

high-order accurate scheme;JFNK method;computational efficiency;scalability

A

：TP301

10.3778/j.issn.1673-9418.1512054

*The National Natural Science Foundation of China under Grant Nos.11272352,61561146395(国家自然科学基金);the Basic Research Program of National University of Defense Technology under Grant No.ZDYYJCYJ20140101(国防科技大学科研计划重大应用基础研究项目);the Open Research Program of State Key Laboratory of Aerodynamics of China under Grant No.SKLA20140104 (空气动力学国家重点实验室开放课题).

Received 2015-11,Accepted 2016-01.

CNKI网络优先出版:2016-01-27,http://www.cnki.net/kcms/detail/11.5602.TP.20160127.1635.004.html