0 引言

随着计算流体力学(CFD)的高速发展，其在航空航天工程中的应用越来越广泛。为了适应航空航天工程对于精细模拟的更高要求，目前CFD已经发展到千万量级的计算网格阶段。非常密的网格给数值计算方法的收敛性带来巨大的挑战，尤其是高雷诺数湍流问题和非定常问题。同时，高性能计算机硬件水平也不断提高，使上万CPU核心的并行计算成为可能。这就需要CFD高效且鲁棒地处理航空航天工程中的复杂外形、适用于最新的高性能计算机及未来的发展^[1-2]。综上所述，提高CFD收敛性和并行计算效率的研究有着重要理论研究意义和工程应用价值。

全隐式方法是CFD提高收敛性和计算效率的重要途径，全隐式方法的并行化是CFD研究难点之一。Luo等^[3]最早将无矩阵LU-SGS方法作为GMRES的预处理算子，在共享内存系统中采用OpenMP实现了算法的并行化。国内阎超等^[4]在混合网格上采用OpenMP实现了GMRES+LU-SGS的并行化。上述的算法均有非常高的收敛效率，但由于共享内存系统的限制，并行效率并不高。Zingg等^[5-6]采用近似Schur并行预处理算子实现了块结构网格上GMRES的分布式并行化，采用无矩阵方法处理矩阵向量积，取得了非常好的收敛效果和并行效果，在一定条件下出现了加速比的超线性，NASA CRM模型千万量级网格湍流模拟仅需3000步收敛。同一算例半隐式的LU-SGS方法需要数万步的迭代才能收敛^[7]。燕振国等^[8]将GMRES应用于高阶耗散紧致格式，张健等^[9]基于PETsc科学计算工具包建立三维混合网格分布式并行GMRES算法，龚小权等^[10]将GMRES应用于间断Galerkin有限元, 均取得了不错的效果。综合分析国内外研究的前沿，在混合网格上GMRES的分布式并行化是隐式方法研究的难点。需要恰当的处理矩阵运算并行化、预处理算子并行化、通讯模型、网格重排序和负载均衡，才能提高GMRES方法并行效率。近年来，计算机硬件高速发展，可用内存越来越大，大规模问题中存储雅可比矩阵成为可能。全隐式方法中计算和存储雅可比矩阵，并采用并行GMRES方法求解线性方程组是提高CFD收敛性的主要途径。

本文采用变子空间数GMRES算法，结合Hybrid LU-SGS方法作为并行预处理算子，引入更鲁棒的湍流模型，结合不同网格排序方法，在MPI并行环境中，发展出一套计算和存储雅可比矩阵的分布式并行全隐式计算方法。在南京航空航天大学自主开发的OVERU软件^[11]中应用该方法，经验证具有非常高的收敛性和并行效率，适用于定常和非定常流问题的数值模拟。

1 控制方程与离散方法

控制方程为任意拉格朗日欧拉(ALE)形式的三维雷诺平均N-S方程(RANS)，湍流模型为一方程SA模型，本文采用的是一种新型的Negative-SA模型^[12]，该模型中${\mathit{\tilde v}}$允许为负值，以抑制网格质量较差时的非物理解，可以改善整体隐式方法的鲁棒性，提高迭代中CFL数的取值，从而间接达到提高全隐式方法收敛性的目的。Negative-SA模型方程为：

当${\mathit{\tilde v}}$≥0时

$ \begin{aligned} \frac{\partial \widetilde{\nu}}{\partial t}+u_{j} \frac{\partial \widetilde{\nu}}{\partial x_{j}}=c_{b 1}\left(1-f_{t 2}\right) \widetilde{S} \widetilde{\nu}-\\ \\\left[c_{w 1} f_{w}-\frac{c_{b 1}}{\kappa^{2}} f_{t 2}\right]\left(\frac{\widetilde{v}}{d}\right)^{2}+\\ \frac{1}{\sigma}\left[\frac{\partial}{\partial x_{j}}\left((\nu+\widetilde{\nu}) \frac{\partial \widetilde{v}}{\partial x_{j}}\right)+c_{b 2} \frac{\partial \widetilde{\nu}}{\partial x_{i}} \frac{\partial \widetilde{\nu}}{\partial x_{i}}\right] \end{aligned} $

(1)

当${\mathit{\tilde v}}$＜0时

$ \begin{aligned} \frac{\partial \widetilde{\nu}}{\partial t} &+u_{j} \frac{\partial \widetilde{\nu}}{\partial x_{j}}=c_{b 1}\left(1-f_{t 2}\right) \mathit{\Omega} \widetilde{\nu}+c_{w 1}\left(\frac{\widetilde{\nu}}{d}\right)^{2}+\\ & \frac{1}{\sigma}\left[\frac{\partial}{\partial x_{j}}\left(\left(\nu+\widetilde{\nu} f_{n}\right) \frac{\partial \widetilde{\nu}}{\partial x_{j}}\right)+c_{b 2} \frac{\partial \widetilde{\nu}}{\partial x_{i}} \frac{\partial \widetilde{\nu}}{\partial x_{i}}\right] \end{aligned} $

(2)

式中各符号与原始SA模型一致。其思想是当${\mathit{\tilde v}}$≥0求解原始SA模型，当出现非物理解${\mathit{\tilde v}}$＜0时求解式(2)，式(2)源项中破坏项c_w1(${\mathit{\tilde v}}$/d)²与原SA模型符号相反，通过迭代求解使得${\mathit{\tilde v}}$重新更新为正值。

OVERU软件采用非结构混合网格有限体积方法求解控制方程。空间离散采用格点格式，可处理多种单元类型的混合网格。基于多种常用的通量差分格式HLLC、ROE、HLLEM和通量分裂格式AUSMPW+等，通过U-MUSCL方法进行高阶重构，限制器采用Venkatakrishna限制器。非定常计算采用双时间步长方法进行二阶精度的时间离散，内迭代采用本文所发展的全隐式方法。

2 并行全隐式方法 2.1 GMRES方法

全隐式方法是定常、非定常CFD计算中提高收敛性和计算效率有效途径之一。其中的Krylov子空间迭代法中的GMRES方法^[13]是一种高效的解线性方程组方法，结合不同预处理算子在CFD数值计算领域广泛使用。本文的分布式并行GMRES方法的架构如图 1所示，含全隐式方法的并行化实现及优化。

在重启型GMRES算法的基础上，本文采用变子空间数GMRES算法，其初始化过程为v₀=R－AΔW₀, r₀:=P^－1v₀, β:=‖r₀‖, v₁:=r₀/β。初始化后，计算流程如下：

$ \begin{array}{l} \begin{array}{*{20}{l}} {01:}&{{\rm{ for }}\;j = 1, m}\\ {02:}&{\;\;\;\;\;{\mathit{\boldsymbol{w}}_j} = {\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}}\Delta \mathit{\boldsymbol{W}}}\\ {03:}&{{\rm{ for }}\;i = 1, j}\\ {04:}&{\;\;\;\;\;{h_{ij}} = \left( {{\mathit{\boldsymbol{w}}_j}, {\mathit{\boldsymbol{v}}_i}} \right)} \end{array}\\ \begin{array}{*{20}{l}} {05:}&{{\mathit{\boldsymbol{w}}_j} = {\mathit{\boldsymbol{w}}_j} - {h_{ij}}{\mathit{\boldsymbol{v}}_i}}\\ {06:}&{{\rm{ end }}}\\ {07:}&{{h_{j + 1, j}} = {{\left\| {{W_j}} \right\|}_2}}\\ {08:}&{\;\;\;{v_{j + 1}} = {\mathit{\boldsymbol{w}}_j}/{h_{j + 1, j}}}\\ {{{09}_:}}&{\;\;\;{z_{1:j}}: = {{\min }_z}{{\left\| {\beta e - {H_{1:i}}{z_{1:j}}} \right\|}_2}}\\ {10:}&{{\rm{ if }}\quad {{\left\| {\beta e - {\mathit{\boldsymbol{H}}_{1:j}}{\mathit{\boldsymbol{z}}_{1:j}}} \right\|}_2} < \varepsilon {\rm{ }}\;{\rm{exit }}\;{\rm{restart }}}\\ {11:}&{{\rm{ end }}}\\ {12:}&{\Delta \mathit{\boldsymbol{W}} = \Delta {\mathit{\boldsymbol{W}}_0} + {\mathit{\boldsymbol{v}}_m}{\mathit{\boldsymbol{z}}_m}}\\ {13:}&{\Delta {\mathit{\boldsymbol{W}}_0} = \Delta \mathit{\boldsymbol{W}}}\\ {14:}&{{\rm{ restart}}\;\;\;{\rm{ GMRES}}\;\;\;\left( {{\rm{goto}}\;\;\;{\rm{01}}} \right){\rm{ }}} \end{array} \end{array} $

(3)

重启型GMRES中通过残差收敛判据确定是否需要GMRES重启，每个重启过程的子空间数不变均为m。通过增加残差极小值问题的求解次数(式(3)中增加第10步)，提前判断是否满足残差收敛判据。当满足收敛判据时，即使子空间数小于m也退出GMRES内循环，从而减少子空间迭代数，提高计算效率。存储空间方面，GMRES方法需要(m+2)×n_nodes×n_eqn×8字节内存空间。其中m为Krylov子空间数，n_nodes为网格点数，n_eqn为方程数，如三维Euler方程为5。

变子空间数GMRES算法，可以在基本不影响收敛性的前提下，减小了每次迭代计算时间，从而提高整体计算效率。变子空间数GMRES算法相比于原重启型GMRES算法在迭代启动若干步的鲁棒性上有一定程度的提高，改善了全隐式方法的启动问题。

2.2 雅可比矩阵计算方法

隐式方法在处理雅可比矩阵时可分为无矩阵方法和计算矩阵方法。其中无矩阵方法一般采用数值差分方法方法直接得到矩阵向量积，不需要计算和存储雅可比矩阵，节约内存空间。有矩阵方法中雅可比矩阵需要提前计算并存储，虽然增加了内存需求和计算量，但是矩阵方法比无矩阵方法具有更好的收敛性。

全隐式方法中雅可比矩阵以双精度实数存储需要(n_nodes+2n_edges)×n_eqn²×8字节内存空间，处理大规模问题时受到内存的限制往往采用无矩阵方法。近几年计算机硬件的发展，内存空间不断增大，价格也不断降低。内存空间已经逐渐不再制约有矩阵方法的应用。雅可比矩阵的计算方法可通过人工编程或自动微分来实现。本文采用人工编程方法实现雅可比矩阵的计算，人工编程方法可以带来较优化的代码和较高的效率。雅可比矩阵由无黏项和黏性项组成，无黏项采用Van-Leer通量分裂雅可比矩阵，黏性项则基于薄层N-S假设。为提高计算效率，实际计算中雅可比矩阵每隔若干次迭代重新计算一次，在这几次迭代中冻结。采用这个策略在不影响收敛性前提下，减少整体计算时间。

无矩阵方法中的矩阵向量积本文采用Luo等^[3]的近似方法，如式(4)所示，简化通量雅可比计算，可节约差分近似中两次残差的计算量:

$ \boldsymbol{A} \Delta \boldsymbol{W} \approx \Delta F=F(\boldsymbol{W}+\Delta \boldsymbol{W})-F(\boldsymbol{W}) $

(4)

2.3 GMRES并行化

基于MPI非阻塞通讯模型，采用分区并行策略实现GMRES的分布式并行化。雅可比矩阵的对角项按点存储，非对角项按边存储，均存储在局部内存中。通过分区边界上的虚拟点进行数据通信，GMRES中矩阵向量积(式(3)中第2步)仅需相邻分区之间的通信，减小通信时间。GMRES中向量内积(式(3)中第4、7步)在各分区计算完成后进行全局归约操作，得到全局Hessenberg矩阵，通信量极小。采用Metis库进行分区实现并行负载均衡，分区后可以保证各分区的网格点数基本一样，同时分区交界面分割的边数尽量少。可扩展性是并行计算方法的重要指标^[14]，非结构网格分区并行实现方法的可扩展性与分区方法的选取有直接关系，问题规模和分区数等比例增大后，如果分区交界面分割的边数基本不变，则并行方法的可扩展性会很好。

Wissink等^[15]在DP-LUR基础上进行改进，发展出并行化的Hybrid LU-SGS方法。该方法的本质是在各处理器的网格分区内部使用LU-SGS方法，进行Gauss-Seidel迭代；在网格分区的边界上，需要数据通信的网格点上使用DP-LUR方法，进行Jacobi迭代。Hybrid LU-SGS方法具有较高的并行效率，且计算时间比DP-LUR方法的计算时间减少约45%。相比其它预处理算子Hybrid LU-SGS效率更高，因此本文选择Hybrid LU-SGS作为GMRES方法中的预处理算子P(式(3)中第2步)，计算流程如下：

交换分区边界数据，然后边界点上有:

$ \boldsymbol{w}_{j}=\boldsymbol{P}^{-1} \boldsymbol{A} \Delta \boldsymbol{W}=\boldsymbol{D}^{-1} \boldsymbol{A} \Delta \boldsymbol{W} $

(5)

分区内部进行LU-SGS向前扫和向后扫:

$ {(\mathit{\boldsymbol{D}} + \mathit{\boldsymbol{L}})\Delta {\mathit{\boldsymbol{W}}^*} = \mathit{\boldsymbol{A}}\Delta \mathit{\boldsymbol{W}}} $

(6)

$ {(\mathit{\boldsymbol{D}} + \mathit{\boldsymbol{U}}){\mathit{\boldsymbol{w}}_j} = \mathit{\boldsymbol{D}}\Delta {\mathit{\boldsymbol{W}}^*}} $

(7)

2.4 网格重排序

本文所发展的数值方法中的网格重排序包含两个方面：

(1) 非结构网格LU-SGS算法中的网格重排序^[3-4]。经过分层和染色两个步骤将网格重排之后，保证每一层中的点与本层其它点不相连、与任一点相连的上下对角点的个数相当，目的是改善上三角矩阵U和下三角矩阵L的平衡性，构造非结构网格超平面，提高LU-SGS方法的效率。本文采用全局排序方法将分层信息发送到各计算分区。

(2) 网格节点序号的重排序^[16]。对网格节点序号进行重排序，其目的是减小稀疏矩阵的带宽，使矩阵中非零项靠近主对角项，减小矩阵存储空间，提高解线性方程组的效率。CFD计算中网格节点序号的重排序还可以减小计算过程中缓存未命中几率，提高计算效率。本文采用RCM(Reverse-Cuthill-McKee)方法进行网格节点序号的重排序。数值试验表明经过RCM排序后可提高约10%的计算效率。

3 收敛特性验证 3.1 二维定常流动

本节将对比不同时间离散方法在二维平面无黏流和NACA0012翼型黏性绕扰流算例中的收敛性。二维平面无黏流计算网格为三角形网格，如图 2，共1426个网格点。图 3给出了残差收敛到机器精度时显式方法(RK)、LU-SGS方法、无矩阵GMRES方法和矩阵GMRES方法(MGMRES)的对比。算例中最大子空间数m=10，GMRES重启次数为2，GMRES内残差收敛判据ε＜0.1，本文后续算例均采用此设置。三种隐式方法均在1000步内收敛到机器精度，收敛性大大优于显示方法，收敛性最好的是MGMRES方法。

NACA0012翼型黏性扰流计算网格采用四边形单元，共23万网格点，如图 4。计算条件为：来流马赫数Ma_∞=0.15，迎角α=10°，雷诺数为Re=6×10⁶。该算例中数值方法的收敛性受到密网格和高雷数严苛的考验。LU-SGS方法CFL数取1×10⁶，GMRES和MGMRES方法均没有启动问题，CFL数可直接取300。图 5分别给出迭代收敛曲线和CPU时间收敛曲线。可见残差收敛6个量级，LU-SGS方法需要约80 000步，GMRES方法约10 000步，MGMRES方法仅需2500步。计算时间上由于MGMRES每迭代步的计算时间较长，迭代初期其计算效率接近或低于GMRES方法。迭代后期MGMRES的收敛速度优势明显，仅需GMRES方法不到二分之一的计算时间。无矩阵GMRES方法在计算效率较LU-SGS方法提高约20%。需要特别说明的是收敛6个量级主要是为了对比不同方法的收敛性，实际计算中并不需要如此大量的迭代。表 1给出了本文计算得到的升阻力系数与文献[17]中不同程序计算结果的对比，符合得很好。

3.2 翼身组合体扰流

对DLR-F6翼身组合体^[18-19]有黏和无黏绕流进行计算分析，其中无黏计算网格采用四面体单元约90万网格点503万单元，黏性网格采用六面体单元，约580万网格点，如图 6。计算条件为：来流马赫数M_∞=0.75，迎角α=0.52°。黏性计算中雷诺数Re=3×10⁶。LU-SGS方法CFL数取1×10⁶，GMRES方法CFL数取300，MGMRES方法无黏算例CFL数取300。MGMRES方法在黏性算例遇到了启动问题，本文种采用变CFL数取值的策略解决启动问题，如下式所示：

$ \mathrm{CFL}=\left\{\begin{array}{ll}{10} & {, \quad n=1} \\ {10 \times 1.2 n-1} & {, \quad 1<n<19} \\ {300} & {, \quad n \geqslant 19}\end{array}\right. $

(8)

首先进行网格节点序号重排序方法的验证，图 7给出了采用网格节点序号的重排序方法(RCM)和不采用网格节点序号的重排序方法(NO-RCM)的对比，其中TET表示无黏扰流算例，HEX为黏性扰流算例。可见RCM方法减小了计算过程中缓存未命中几率，经过RCM排序后可提高约10%的计算效率。

图 8给出了黏性扰流算例中重启型GMRES方法和本文改进的变子空间数GMRES方法残差收敛曲线的对比，两者收敛效率一致。为达到GMRES内残差收敛判据，当重启型GMRES方法中一般需要1次重启，共10+10次子空间运算，而变子空间数GMRES方法在重启后子空间数会减少，共10+n次子空间运算(n≤10)。变子空间数GMRES方法较重启型GMRES方法，计算时间减少约15%。

图 9给出了三种隐式方法在进行黏性扰流计算时的迭代收敛曲线和CPU时间收敛曲线。其中MGMRES方法的收敛性最好，仅需800步左右残差下降3个量级。残差下降3个量级MGMRES方法较GMRES方法计算时间减少约37%，整体计算时间是LU-SGS方法的四分之一不到。并且MGMRES方法在残差下降3个量级后收敛性依然较好，而GMRES方法收敛曲线出现一定程度波动。图 10给出了升力系数的迭代和CPU时间收敛曲线。从升力系数上判断MGMRES方法收敛效率也是最高，略高于GMRES方法，计算效率约为LU-SGS方法的4倍多。

针对第六届AIAA阻力测试会议中的NASA CRM翼身组合体构型，采用NASA提供的混合网格约3000万网格点1.2亿网格单元，验证本文的方法在数千万量级网格上的收敛性。计算条件为：来流马赫数Ma_∞=0.85，迎角α=2.75°，雷诺数Re=5×10⁶。图 11给出了三种隐式方法在进行黏性扰流计算时的升力系数收敛曲线。其中MGMRES方法升力系数收敛仅需1000步左右，而GMRES方法需要约1600步。相同计算状态及网格条件下，文献[20]中NASA的FUN3D软件采用全隐式方法需要约3000步收敛。本文的方法在千万量级网格上的也具有很好的收敛性。

3.3 战斗机非定常气动弹性模拟

本节采用F16战斗机外形的非定常气动弹性数值模拟算例，考察所发展的隐式方法在三维非定常黏性扰流计算中的效果。计算网格采用混和网格，图 12为机身物面网格和对称面网格，约214万网格点。计算条件为：来流马赫数Ma_∞=0.8，迎角α=4.0°，CFL数取值与3.2节一致。结构有限元模型基于ECERTA项目^[21]提供的结构模态数据。在时域内耦合求解非定常N-S方程和结构模态方程，实现飞机气动弹性的数值模拟。其中流固界面插值方法为薄平板样条法(TPS)，动网格方法采用高效的Delaunay图映射方法。双时间步长方法中内迭代残差下降3个量级为判断其收敛的依据。

图 13给出了三组双时间步长内迭代残差曲线，MGMRES方法仅需约7步迭代残差下降3个量级，GMRES方法需要约8步，LU-SGS方法则需要约20步。收敛性仍然是MGMRES方法最好。计算时间上GMRES方法略快于MGMRES方法，较LU-SGS方法减少约23%。需要注意的是，此时间仅统计了流场计算的时间，未考虑气动弹性计算中网格变形、流固界面插值等时间。

图 14给出了采用MGMRES方法计算得到的涡量等值面及压强云图，符合扰流基本规律。图 15给出了前四阶模态广义位移随时间的变化曲线，各阶广义位移振幅逐渐减小，在此状态飞机未发生颤振。

3.4 旋翼前飞流场

将本文方法发展到非结构重叠网格上，通过旋翼非定常流模拟研究方法的收敛特性。采用2.5m级单旋翼模型，包括4片桨叶。旋翼实度0.077，翼型采用NACA0012，弦长0.072 m，几何扭转0°，前飞马赫数0.127。计算网格方面，如图 16所示采用的非结构重叠网格系统，共包含6个子网格，其中两块为背景网格点数分别为601万和139万，每片旋翼叶片的网格随旋翼运动，网格点数为107万，总网格数约1200万。

双时间步长方法中内迭代残差下降4个量级为判断其收敛的依据。在内迭代中MGMRES方法遇到了鲁棒性问题，非定常CFL最大仅可取到80，而LU-SGS方法可取1×10⁶。图 17给出了采用MGMRES方法计算得到的Q等值面，其中桨涡干扰现象明显。图 18给出了两组内迭代残差曲线，MGMRES方法仅需约35步迭代残差下降4个量级，LU-SGS方法则需要约107步。收敛性仍然是MGMRES方法最优。计算时间上MGMRES方法LU-SGS方法减少约40%。但是，内迭代残差下降3个量级可满足一般的非定常计算需要，此时LU-SGS方法在计算时间上还要少于MGMRES方法。主要原因是此算例中MGMRES方法的CFL数较小，当达到GMRES内残差收敛判据ε＜0.1时，子空间数仅有1~3个，降低了GMRES方法的收敛速度。因此提高GMRES方法的鲁棒性的研究在未来仍然是很有必要的。

4 并行特性验证

本节采用3.2节的翼身组合体绕流算例进行本文发展方法的并行特性验证。图 19为各隐式方法内存空间需求，无矩阵GMRES方法比LU-SGS方法内存需求增加有限。MGMRES方法需要的内存最多，黏性算例在580万网格上需要约13GB内存。MGMRES方法平均每千万网格需要约22GB内存空间。在目前硬件条件下，采用MGMRES方法单台工作站就能够满足千万量级网格的内存需求，高性能分布式集群则可进行更大规模网格的计算。

在国家超级计算天津中心的天河1A高性能计算机上，对所发展方法的并行特性展开研究。图 20为采用MGMRES方法、1和240个计算核心(NP)计算得到的机翼上四个站位(η)压强系数分布与实验值的对比，两者符合得较好，验证了并行方法的正确性。

图 21给出了Hybrid LU-SGS方法和LU-SGS方法残差收敛曲线。Hybrid LU-SGS方法在分区边界网格点上使用DP-LUR方法而LU-SGS方法对分区边界不做处理(图中w/o表示)。可见在分区较少时两种方法残差收敛曲线非常接近。当240个分区时，分区边界网格点不做处理的并行方法在4000步左右计算发散。因此，在进行较多分区并行计算时，采用Hybrid LU-SGS方法是必要的。图 22给出了Hybrid LU-SGS预处理算子和MGMRES方法并行加速比。两种方法的加速比非常接近理想值，加速效率非常高，MGMRES方法的加速比略低于Hybrid LU-SGS方法。特别的，在240个计算核心的并行计算中均出现了加速比超线性(并行效率超过100%)。一般认为，在负载均衡和极端通信时间条件下，加速比超线性主要是因为计算机CPU缓存在更多分区时利用率更高。

5 结论

本文基于变子空间数GMRES算法，结合Hybrid LU-SGS方法作为并行预处理算子，引入更鲁棒的湍流模型，结合不同网格排序方法，发展出一套计算和存储雅可比矩阵的分布式并行全隐式计算方法。利用南京航空航天大学自主开发的OVERU软件，在天河1A高性能计算机上对若干二维和三维扰流问题进行了数值模拟，经验证具有非常高的收敛性和并行效率，结论如下：

(1) 变子空间GMRES方法与重启GMRES方法具有相同的收敛特性，计算时间可较原方法减少约15%。计算和存储雅可比矩阵的全隐式方法收敛性优于无矩阵方法，较LU-SGS方法提高20%~200%计算效率。采用RCM网格重排序方法优化缓存命中可将计算效率再提高10%。

(2) 矩阵GMRES方法每千万网格约需要22GB的内存空间，当前硬件水平下可以满足大规模网格计算的内存需求。

(3) 分布式并行GMRES方法结合Hybrid LU-SGS预处理算子，所发展的全隐式方法并行效率非常高。特别的，在240个计算核心出现了加速比的超线性。

本文研究仅限于CPU并行计算，后续研究应面向未来E级超算的CPU+GPU(或加速卡)的异构架构，开展异构架构下全隐式并行算法研究。