• 返回主页
  • CFD中的多重网格求解器

    2016.05.19,Update:增加C++求解代码;
    众所周知,在对控制方程进行离散后,会产生一个离散的数学方程组,这个方程组通常为稀疏的。其大小和复杂性跟维度、方程个数、离散格式有关。虽然我们可以使用任意一个方程组求解技术来求解这个方程组,但是目前的限制主要为计算机能力。目前的求解方法主要有直接求解法和迭代求解法。直接求解典型的方法为高斯消去法。在使用高斯消去法的时候,需要对求解系统进行$N*N*N$次操作($N$为未知数数量),并且需要利用内存存储$N*N$个系数。迭代方法即为通过一系列迭代的方法获得最后的解。典型的迭代方法为Jacobi和Gauss-Seidel点迭代法。在迭代求解中,一般在每次迭代中操作$N$次即可,但是总的迭代数是不确定的。

    目前在大部分CFD软件中,均采用迭代求解法来求解方程离散后的稀疏线性系统。然而,Jacobi和Gauss-Seidel点迭代法的思想虽然很简单,但是在处理大型方程组的时候收敛速率是非常慢的。其可以较快的消除短波误差分量,但是对长波误差分量的消除非常缓慢,因此他们并不适用于CFD求解。得益于多重网格求解技术,Gauss-Seidel点迭代法目前也在商业CFD软件中广泛应用。

    多重网格方法的求解思想为首先在较细的网格上进行迭代来消除短波误差分量,然后再在粗网格上迭代消除长波误差分量。通过这种在不同粗细的网格上消除不同波长的误差分量的方法,可以快速的把各种波长的分量消除。在每个不同粗细量级的网格上,可以使用任何一种迭代求解技巧。

    本文从一个包含三个未知数的方程组开始,简单介绍Gauss-Seidel点迭代法,然后介绍多重网格求解思想的基本步骤。最后使用基于几何的多重网格方法(Geometric multigrid, GMG)求解一个包含20个未知数的方程组并和传统的Gauss-Seidel点迭代法进行收敛速度对比。

    一、Gauss-Seidel点迭代法

    Gauss-Seidel点迭代法在大部分矩阵、数值分析教科书中均有提及。在这里用一个说明性的例子阐述Gauss-Seidel点迭代法的基本思想。

    假设我们有如下的方程组: \begin{matrix}\tag{1} 2x_1+x_2+x_3=7 \\ -x_1+3x_2-x_3=2 \\ x_1-x_2+2x_3=5 \end{matrix} 对方程移项: \begin{matrix}\tag{2} x_1=(7-x_2-x_3)/2 \\ x_2=(2+x_1+x_3)/3 \\ x_3=(5-x_1+x_2)/2 \end{matrix} 假定初值为: \begin{matrix}\tag{3} x_1^0=0 \\ x_2^0=0 \\ x_3^0=0 \end{matrix} 将方程(3)中的$x_1^0$代入到方程(2)中有: \begin{matrix}\tag{4} x_1^1=(7-x_2^0-x_3^0)/2=(7-0-0)/2=3.5 \end{matrix} 将方程(4)求解的$x_1^1$和方程(3)中的$x_2^0$代入到方程(2)中有: \begin{matrix}\tag{5} x_2^1=(2+x_1^1+x_3^0)/3=(2+3.5+0)/3=1.8333 \end{matrix} 将方程(5)求解的$x_2^1$和方程(4)求解的$x_1^1$代入到方程(2)中有: \begin{matrix}\tag{6} x_3^1=(5-x_1^1+x_2^1)/2=(5-3.5+1.8333)=1.6667 \end{matrix} 这样,我们有Gauss-Seidel点迭代法求解的第一次初始值: \begin{matrix}\tag{7} x_1^1=3.5 \\ x_2^1=1.8333 \\ x_3^1=1.6667 \end{matrix} 在接下来的迭代过程中,$x^1$的值即可当做初始值进行下一次迭代。针对方程(1),10次迭代业已收敛。

    二、多重网格求解思想

    多重网格求解技术需要在不同粗细量级的网格上对同一个原始的PDE方程进行求解。考虑最细的原始网格均匀分布且节点距离为$h$,假设粗网格节点距离为$2h$,$4h$,$8h$...。首先我们需要在间距为$h$的网格上进行若干次迭代,然后跳跃至$2h$的网格上进行下一轮迭代。如图所示,最终在$h$网格上获得一个加速收敛的解。其具体步骤为:
    多重网格
    图1. 多重网格V-cycle示意图,其中$h$为网格节点距离,横线表示不同粗细量级的网格,圆圈内的数字表示每层网格上的迭代数

    1. Start multigrid iterations:

    \begin{matrix}\tag{8} A_h x = b \end{matrix} 其中$A_h$表示在间距为$h$的网格上对原始PDE离散后的矩阵系数。此处可以使用任意一种迭代技术求解,如使用Gauss-Seidel点迭代法迭代10步。在此步,可以快速的消除短波误差分量。迭代后我们有解$y_h$,以及残差$r_h$。

    2. Restriction:

    2.1:基于间距为$h$的网格上的残差$r_h$,需要通过Restriction对间距为$2h$的网格上的残差进行操作。一种简单的做法即为平均插值方法或带权插值方法(非均分网格)。需要注意的是,$r_h$共有$N$个分量,$r_{2h}$则有$2N$个分量。通过Restriction操作后,假设我们有3套粗网格,需要执行三次Restriction,最后有:$r_{2h}$,$r_{4h}$,$r_{8h}$。

    2.2:在前一步,我们通过在间距为$h$的网格上对原始PDE离散后得到系数矩阵$A_h$,但是并不知道$A_{2h}$,$A_{4h}$,$A_{8h}$的具体值。多重网格可以使用基于几何(Geometric)或者基于代数(Algebraic)的方法来求出$A_{2h}$等矩阵系数。这俩种方法也即GMG和AMG。

    2.3:在获得$r_{2h}$,$r_{4h}$,$r_{8h}$后,我们可以通过$A_{2h}$,$A_{4h}$,$A_{8h}$对误差$e_{2h}$,$e_{4h}$,$e_{8h}$进行迭代求解。通常并不需要很多步骤,简单几次即可。求解后我们有最终的$e_{2h}$,$e_{4h}$,$e_{8h}$。

    2.4:通过误差$e_{2h}$,$e_{4h}$,$e_{8h}$对残差$r_{2h}$,$r_{4h}$,$r_{8h}$进行更新有最后的$r_{2h}'$,$r_{4h}'$,$r_{8h}'$。

    3. Prolongation:

    3.1:基于$e_{8h}$,插值计算$e_{4h}'$。然后对$e_{4h}$和$e_{4h}'$加和我们有新的$e_{4h}''$。以$e_{4h}''$作为初值,使用任意的迭代技术对$A_{4h} e_{4h}''=r_{4h}'$进行光顺(smooth)求解后有光顺的误差$e_{4h}'''$。在下文中用$e_{4h}$代替$e_{4h}'''$。

    3.2:同理,基于$e_{4h}$(也即上文的$e_{4h}'''$),重复步骤3.1,我们有光顺后的误差$e_{2h}'''$。在下文中用$e_{2h}$代替$e_{2h}'''$。

    3.3:同理,基于$e_{2h}$,最终我们有$e_{h}'''$。

    4. Correction:

    计算最终的解:$y=y+e_{h}'''$。至此我们完成了多重网格的一次sweep。重新进行第1步。

    需要提及的是,在Prolongation的过程中,通常对误差进行光顺(smooth),Gauss-Seidel是一个非常好的光顺器。下面我们使用基于几何的多重网格(GMG)方法求解一个包含20个未知数的方程组来详细解释多重网格方法的具体实现步骤。

    代码下载:使用C++对一个简单的算例进行多重网格求解,直接编译即可。 点击下载

    代码解析待续


    东岳流体
    info@dyfluid.com