LU-SGS求解器





  • 或许你可以试试看 :mianmo:



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



  • @赵一铭
    最近看了一下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,不过也分点隐、线隐之类的玩意儿,看上去整得挺高级。



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

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

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

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

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



  • @李东岳
    @李东岳

    可以,但是精度个人感觉取决于对角占优的程度。原方程是:
    $$(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实现来看需要对矩阵系数进行操作,里面数据依赖关系决定了可能还是得把矩阵组装出来。



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

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



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

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



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

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



  • @程迪

    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是一个非迭代的求解系统是正确的,问题是从我测试的那个矩阵求解来看,结果只能说是凑合。还是说我的理解有偏差?



  • @李东岳

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

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



  • @李东岳

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

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

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

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



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

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

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



  • @李东岳

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

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



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



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



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



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



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

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



  • @程迪 怎么jacobian-free,在JVNK中也是以matrix-free实现的。AIAA很多论文都是如此表述的。matrix-free的LU-SGS还真心free的飞起



  • @youmengtian
    Jacobian free是不算Jacobian,JFNK的方法压根就没出现Jacobian,既不存,也不算;Matrix free主要是不存系数矩阵,至于这个矩阵是不是Jacobian无所谓的,但是一般多少都有点儿关系。Matrix free不代表不算矩阵元素,preconditioner免不了要算全部或者部分矩阵元素的,但是往往算完之后可以扔掉一些,不用都存下来。SGS就是一例,按照LU-SGS的搞法,就是Matrix free的,计算过程中用到了某些系数,就算一下,用完了就扔了,下次要用再算。



  • @程迪 是这么说,目前的问题就是还是对于编写这类耦合隐式算法没有思路。OpenFOAM也没有相关可以借鉴的代码。不知您在这方面是否已经有相关工作开展了?可以进行进一步交流。



  • @youmengtian
    没,网上有个那个Chun Shen的代码,AETK,我正在看。
    https://github.com/chengdi123000/AETKv1



  • 不好意思,最近除了上班,就是在准备课程资料和写CFD界,一直没空看这个LUSGS。等我忙完了我看看

    • 11楼我的测试初步感觉LUSGS一步求解稀疏矩阵求解精度不高,大约残差在0.01左右,这在SIMPLE算法里面是可以容忍的,因此本身SIMPLE就是迭代求解的过程。但是对于非迭代求解的PISO算法,瞬态问题需要每个时间步都收敛,在这种情况下,LUSGS如何能保证准确性?

    这是目前我需要理解的地方,谁敢兴趣也可以交流一下。
    LUSGS比较有意思,等搞明白了发个CFD界,毕竟之前都是迭代求解矩阵,LUSGS一步求解,我觉得不错。



  • @李东岳 感兴趣呀😄



  • @李东岳

    lusgs主要算空动问题,和piso算法可能从来没有过交集吧。

    而且计算代价比较低可以当预条件用。


  • OpenFOAM讲师

    @程迪 请问什么叫空动问题?



  • 估计是空气动力学问题。如果是密度基求解器的话,每个时间步一次性收敛是更重要的了。



  • @李东岳
    我觉得时间推进是非线性问题,咋可能一步收敛。只要
    1.对时间精确的求解来说,不影响所需要的求解精度;
    2.对不要求时间精确的,只要不影响下一次迭代就行;



  • @yhdthu aerodynamics



  • LUSGS的好处是可以使用更大的时间步长,但是显式计算的RK4也可以做到CFL~0.8。在二阶程序里使用LUSGS, CFL还是不应该大于1,因此LUSGS失去了它的优势。如果希望计算稳态,LTS也是一种相当有效率的稳态算法。如果LUSGS计算稳态收敛的速度比LTS更快,那么还是值得使用它。



  • 'Implementation of density-based solver for all speeds in the framework of OpenFOAM' 里面介绍的双时间推进看起来不错,虽然计算量大一点,但是应该可以更好的处理流场马赫数跨度很大的情况 (跟压力基础求解器相比)。我一同事刚好在算这种流场,也许可以尝试你写的那种算法。



  • 这里讨论的实在是高深。虽然不懂 LU-SGS 是啥,但最近偶然看到有人分享这样的代码:
    网站
    GitLab



  • @程迪
    您好,有个问题想请教您。看到您在《OpenFOAM的残差定义》一文中提到了残差定义方式:https://chengdi123000.github.io/2018/01/04/OpenFOAM的残差定义/
    804487c1-8422-49c6-a9b1-116b1949a77f-image.png
    关于分母normFactor,有些疑问,不知道为什么要这样定义,有什么数学上的原理吗?
    谢谢!


Log in to reply
 

CFD中文网 2016 - 2020 | 京ICP备15017992号-2