LU-SGS求解器




  • 网格教授 OpenFOAM教授 管理员

    或许你可以试试看 :mianmo:


  • OpenFOAM讲师

    LU-SGS 听名字像线性代数求解器,但是感觉看公式不像是线性代数求解器,反倒是像某种格式。


  • OpenFOAM讲师

    @赵一铭
    最近看了一下LU-SGS的相关资料。我的感想如下:

    LUSGS从算法来讲是近似求解线性代数方程的方法。大概的关系是:
    对于$(L+D+U)x=b$的方程,精确求解可以用Gauss-Seidel (GS)迭代:

    • Gauss-Seidel迭代是:$(L+D)x^{n+1}=b-Ux^n$
      进一步可以发展成Symmetric Gauss-Seidel 对称GS迭代:
    1. $(L+D)x^{*}=b-Ux^n$
    2. $(U+D)x^{n+1}=b-Lx^*$

    这两者是迭代方法,也就是每次迭代残差减小,但是得反复求解多次才能得到解。不过优势是得到的解是线性方程的精确解(在一定残差意义下)。

    对于CFD而言,这样的求解方式有几个问题:

    • 由于其中的数据依赖关系,需要组装LDU矩阵;
    • CFD有非线性误差,所以线性方程解得太准没意义;

    所以有人搞出了LU-SGS。
    从算法来看,其实就是$x^{0}=0$的SGS,但是只迭代一次SGS,得到$x^1$就收工。
    好处是:

    • 由于根据其中的数据依赖关系,不用组装LDU矩阵;
    • 得到的解是不精确的,但是对于对角占优的矩阵,也差不太多了;

    具体使用的时候还涉及另一个技巧:Newton线化法。

    非线性问题线化有两种方法:Newton迭代和Picard迭代。

    • Picard迭代就是OpenFOAM中常用的。$U\cdot U$项被线化成$U^{n} \cdot U^{n+1}$,理论上似乎是一阶收敛的,好处是简单,不用求导。

    • Newton迭代是要求导数的,是$U \cdot U$线化成$U^{n} \cdot U^{n}+2U^{n}\cdot \Delta U^{n}+(\Delta U^{n})^2$,其中$\Delta U^{n} = U^{n+1}-U^{n}$表示增量。然后你就能看到LU-SGS理论里整出的方程是Residual Form的。

      • newton法的好处是收敛快,理论2阶;配合Residual Form的好处是数值误差可以小一些。
      • 额外的好处是:newton法可以玩不精确的导数项(多维时候是Jacobian项),实际LU-SGS也是这么干的。配合Residual Form可以玩出以下花样:
      • 方程右侧的Residual Form反正是显式的,可以用准确一点儿的高精度格式来算;左侧的导数项可以用近似导数(减小计算量,减少数据耦合),只要导数符号上没有大问题,求个稳态解是没啥问题的,收敛速度也可以达到近似二阶。而且收敛到最后稳态解残差是为0的,所以相当于收敛到了一个高精度格式的解。
      • 最后,有个稳态求解器,用那个啥伪时间步可以搞出瞬态的求解器哦。

    具体那文章附送代码,不过好像挺复杂的,没有用coupled solver,直接手工搞的LU-SGS。

    现在高级一点儿LU-SGS似乎是Block LU-SGS,不过也分点隐、线隐之类的玩意儿,看上去整得挺高级。


  • 网格教授 OpenFOAM教授 管理员

    LUSGS从算法来讲是近似求解线性代数方程的方法。

    所以说LUSGS可以用于任何的矩阵求解?比如icoFoam里面的U或者p或者epsilon,也可以用LUSGS求解?

    从算法来看,其实就是x0=0x0=0的SGS,但是只迭代一次SGS,得到x1x1就收工。

    前面的有关Gauss-Seidel (GS)迭代部分很清楚,这句话有点没跟上思路,哈。然后后来的Newton迭代和Picard迭代没太跟上他们和LUSGS的关系,

    楼上好人,感谢分享 :cheeky:


  • OpenFOAM讲师

    @李东岳
    可以,但是精度个人感觉取决于对角占优的程度。原方程是:
    $$(L+D+U)\cdot x = b$$
    这等价于
    $$(L+D)D^{-1}(D+U)\cdot x = b-LD^{-1}U\cdot x$$
    然后LU-SGS求解的是:
    $$(L+D)D^{-1}(D+U)\cdot x = b$$
    LU-SGS与精确线性方程的偏差项是$LD^{-1}U\cdot x$。如果对角占优,$D^{-1}$貌似可以保证这一项比较小,我数值试验发现额外的好处是$LD^{-1}U$是和M有相同的稀疏结构。

    然后你把LU-SGS的公式和SGS$x^0=0$的第一次迭代相比较,可以发现它只是等效于一次$x^0=0$的对称GS迭代而已。

    所以如果你组装矩阵,然后用LU-SGS,还不如直接用SGS多迭代几次,线性方程有时候还是收敛得精确一点儿好。

    而其实LU-SGS最大的好处在于不用组装矩阵。它的数据依赖关系决定了只需要做几次sweep就行了。和OpenFOAM中实现GS的方式差不多。我个人感觉差别在于对于耦合的情况也可以这么干。

    Newton和Picard是两种不同的线性化方法:

    • Picard线化得到的是全量方程,不用求Jacobian,理论一阶收敛;
    • Newton线化得到的是增量方程,具体就是一个Residual Form的方程,至少得给个近似Jacobian,比如用JFNK之类的;

    OF里面的解我看到的都是Picard迭代进行线性化,理论收敛速度慢于Newton迭代线性化。

    我的理解是:LU-SGS相当于给了一个近似的,可以快速求解的Jacobian(两次sweep,不用组装矩阵)。所以必须在采用Newton线化方式的时候才能发挥最大效力。不然没有太多用处。


  • OpenFOAM讲师

    @李东岳
    @李东岳

    可以,但是精度个人感觉取决于对角占优的程度。原方程是:
    $$(L+D+U)\cdot x = b$$
    这等价于
    $$(L+D)D^{-1}(D+U)\cdot x = b-LD^{-1}U\cdot x$$
    然后LU-SGS求解的是:
    $$(L+D)D^{-1}(D+U)\cdot x = b$$
    LU-SGS与精确线性方程的偏差项是$LD^{-1}U\cdot x$。如果对角占优,$D^{-1}$貌似可以保证这一项比较小,我数值试验发现额外的好处是$LD^{-1}U$似乎是和$M$有相同的稀疏结构。

    然后你把LU-SGS的公式和SGS取$x^0=0$的第一次迭代相比较,可以发现它只是等效于一次$x^0=0$的对称GS迭代而已。

    所以如果你组装矩阵,然后用LU-SGS,还不如直接用SGS多迭代几次,线性方程有时候还是收敛得精确一点儿好。

    而其实LU-SGS最大的好处在于不用组装矩阵。它的数据依赖关系决定了只需要做几次sweep就行了。和OpenFOAM中实现GS的方式差不多。所以我觉得对于标量方程可能并没有什么特别的好处。但是对于矢量或者耦合方程而言,组装矩阵可能代价就比较大了。

    举个例子:1M自由度(Degree of freedom, DoF,可以理解为网格数),1个double占8字节=8B,一个标量方程的自变量(比如$p$)占8MB内存,而可压RANS有$\rho, u ,v, w, e, k, eps$ 7个变量,一个时刻的场变量大概是56MB内存,这个是免不了的。这个随自由度线性增长,没有办法的事情。分离求解可以想办法少一点儿。

    但对于方程系数矩阵而言,

    • 如果不用稀疏矩阵,每个标量方程的系数就是1M*1M*8B=8TB内存,鬼都受不了。
    • 用稀疏矩阵,分离求解,比如LDU形式,每个四面体单元算4个面吧,每个面被2个单元使用,所以差不多有$size(L)=size(U)\approx size(D)*2$,大概用(2*2+1)*1M*8B=40MB就存下了。而且随网格数线性增长,分离求解,每次求完一个变量就可以把系数矩阵LDU清空给下一步求解用,所以总的内存消耗量还挺少的;
    • 用稀疏矩阵,耦合求解,LDU每个元素都是小的7*7的block,(按OpenFOAM的尿性通常是湍流单独求解,把其他变量耦合起来也得是5*5的block),每次求解的系数矩阵就是(2*2+1)*5*5*1M*8B=1GB内存。大内存的读、写、元素的计算啥的都拖慢速度。这点对于GPU,众核之类的平台可能更加重要,PC内存现在跟白菜似的。
    • 稀疏矩阵,LU-SGS,每次只要两此sweep遍历所有单元,没有系数矩阵,每次只要大概(2*2+1)*5*5*8B=1MB内存就够了。对角占优时精度也不赖。

    总结:LU-SGS省内存,对于隐式耦合求解来说最合适不过了。

    Newton和Picard是两种不同的线性化方法:

    • Picard线化得到的是全量方程,不用求Jacobian,理论上一阶收敛;
    • Newton线化得到的是增量方程,具体就是一个Residual Form的方程,理论上二阶收敛,但至少得给个近似Jacobian,比如用JFNK之类的;

    OF里面的解我看到的都是Picard迭代进行线性化,理论收敛速度慢于Newton迭代线性化。

    我的理解是:LU-SGS相当于给了一个近似的,可以快速求解的Jacobian(两次sweep,不用组装矩阵,快,省内存)。LU-SGS等价于单步$x^0=0$的SGS这一点也和Newton线化得到的增量形式方程也很契合,和非精确Newton法只需要近似Jacobian更加契合。所以必须在采用Newton线化方式的时候才能发挥最大效力。不然没有太多用处。

    替代品Alternative: Krylov系列求解器也是个省内存狂魔,理论上idr(1)。但是其收敛速度和收敛性取决于矩阵特征值的集中程度,越集中越快,所以似乎又需要预条件(预条件是降低条件数$cond(M)=|\frac{\sigma_{max}(M)}{\sigma_{min}(M)}|​$,其实也差不多就是增加特征值的集中程度)。预条件我不懂,就OF实现来看需要对矩阵系数进行操作,里面数据依赖关系决定了可能还是得把矩阵组装出来。


  • 网格教授 OpenFOAM教授 管理员

    非常非常感谢你的知识分享。内容很大,我得消化消化,暂时我先上班,北京时间周日凌晨我在这里发表一写学后感 :expressionless:

    这种知识和作者本身乐于分享并促进交流,应该有更高的曝光率,我在琢磨以一种什么样的方式,或许可以在我那个网站添加个链接过来,目前东岳流体月PV有20000。


  • OpenFOAM讲师

    @李东岳
    我咋没法编辑答案了。有些地方我想最好修改一下。

    p.s.我觉得有时候最难理解的不是怎么做,而是为什么这么做。



  • 我咋没法编辑答案了。有些地方我想最好修改一下。

    之前有人问问题之后,获得到答案之后就把问题删除了,后来我们就说防止编辑,不准确的可以继续回帖。
    你这个答案我也学习学习… 目前双曲PDE这面对流格式的精准度也是个问题,不知道LU-SGS和我这个有没有关系,当时我只是有这么个想法.


  • 网格教授 OpenFOAM教授 管理员

    @程迪

    SGS有:
    \begin{equation}
    (L+D)x^{ * }=b-Ux^n , (U+D)x^{n+1}=b-Lx^{ * }
    \end{equation}
    给定一个初始$x^0$,可以依次推进。从你提供的LU-SGS公式来看,SGS实际求解的是:
    \begin{equation}
    (L+D)D^{-1}(D+U)\cdot x = b-LD^{-1}U\cdot x
    \end{equation}
    然后忽略掉$-LD^{-1}U\cdot x$,方程变成为:
    \begin{equation}
    (L+D)D^{-1}(D+U)\cdot x = b
    \end{equation}
    这个就是LU-SGS最后求解的系统,不存在一个迭代的过程?是一个非迭代的求解器?

    $-LD^{-1}U\cdot x$这一项忽略掉了是LU-SGS的一个假设,也就是说LU-SGS中认为在对角占优的情况下$-LD^{-1}U\cdot x$足够小。我测试了一下,确实相比原矩阵要小。然后我尝试计算了以下矩阵,结果看起来也还凑合
    0_1497165115792_草图6.png
    上图为通过逆的方法直接求解

    0_1497165122459_草图4.png
    上图为LU-SGS非迭代求解

    咱们先从LU-SGS求解器开始,暂且不讨论效率和线性化,各个击破:sunglasses:
    如果我的理解LU-SGS是一个非迭代的求解系统是正确的,问题是从我测试的那个矩阵求解来看,结果只能说是凑合。还是说我的理解有偏差?


  • OpenFOAM讲师

    @李东岳

    没错,我的理解也是,LU-SGS是一个非迭代的求解器。只有sweep,没有iteration, convergence check之类的东西。

    我尝试过随机生成不对称对角占优的稀疏矩阵,精确解和LU-SGS解的相关性很高。估计是有理论证明的。


  • OpenFOAM讲师

    @李东岳

    这涉及到非线性误差和线性代数求解器误差的相对关系。

    迭代开始前的初始残差基本可以认为是非线性误差,一般线性求解器的收敛准则都设置为下降2个量级,也就是relTol 0.01,甚至有的是relTol 0.1。大部分情况LU-SGS应该是足够了的。

    其实线性代数问题只要愿意迭代,总有收敛的时候,精度只受条件数和机器精度限制,和你采用什么数值方法没关系,只是收敛效率有所区别。但是非线性问题的精度和采用的格式、步长啥的关系都比较大了。

    但是即使线性代数问题解到1e-16,下一个时间步或者非线性迭代步的初始误差又回去了,通常是1e-3左右,和CFL数一般呈正相关,所以线性代数问题解到1e-16没有意义,只要解到1e-5左右就好了,CFL数越大可以考虑继续减少线性问题的求解精度。


  • 网格教授 OpenFOAM教授 管理员

    刚才看了一下一楼那个论文,从方程12,13来看应该确认了LU-SGS是一个非迭代的求解器。测试来看LU-SGS应该可以达到稳态外循环每一步收敛的要求(relTol 0.01或者0.1),但是瞬态的情况下,是需要收敛的。且在可压缩高速密度基求解器中,连外循环都没有。不过文章中说他也可以适用于瞬态算法,看起来LUSGS主要用于高速可压缩流。

    另外,rhoCentralFoam中,以及一些文献中,高马赫数都要采取通量分裂处理,这个原因是什么?
    在我的多相流算法中,也采用了通量分裂处理,不过我们那个的原因只是为了数值稳定,再详细的我就不在这说了有点跑题。可压缩流中不进行通量分裂会如何?

    我再仔细看看那个文章再继续讨论。


  • OpenFOAM讲师

    @李东岳

    高速采用通量分裂是为了在某个特征方向上搞成上风。不然激波、接触间断啥的搞不定。

    不分裂,纯上风也没问题呀,


  • 网格教授 OpenFOAM教授 管理员

    抱歉最近在准备一些课程资料,这个我看看这周末有没有时间详细看看,嘿嘿。



  • 最近着手做matrix-free的LU-SGS实现,本质上是不需要matrix数据结构的。但同时也遇到多重网格如何与之协同工作的问题。目前还只是尝试阶段,希望顺利


  • OpenFOAM讲师

    @youmengtian
    感觉MG和LUSGS好像不是一路的吧。不过LU-SGS当个smoother应该问题不大。



  • @程迪 目前fluent的做法是使用point gs法或ILU(0)做smoother。目前的问题就在无矩阵,或更科学的叫jacobian-free,其本质不再是传统的Ax=b的形式,很难直接用相应的solver。


  • OpenFOAM讲师

    @youmengtian
    明显Jacobian Free和matrix free不是一回事儿吧。JFNK虽然是Jacobian free,但是Krylov没有preconditioner不是很好使,加preconditioner又很难避免要存、算矩阵,至少算矩阵系数是免不了的,不过这个preconditioner反正是approximation,似乎不需要非常consistent。

    而且preconditioner整的不好就是并行的性能瓶颈。