LU-SGS求解器

赵一铭

这个算法有人实现过么？

Implementation of density-based implicit LU-SGS solver in the framework of OpenFOAM

李东岳

或许你可以试试看

程迪

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迭代：
$(L+D)x^*=b-Ux^n $
$(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中，以及一些文献中，高马赫数都要采取通量分裂处理，这个原因是什么？
在我的多相流算法中，也采用了通量分裂处理，不过我们那个的原因只是为了数值稳定，再详细的我就不在这说了有点跑题。可压缩流中不进行通量分裂会如何？

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

程迪

@李东岳

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

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

李东岳

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

youmengtian

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

程迪

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

youmengtian

@程迪目前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整的不好就是并行的性能瓶颈。

youmengtian

@程迪怎么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的，计算过程中用到了某些系数，就算一下，用完了就扔了，下次要用再算。

youmengtian

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

程迪

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

李东岳

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

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

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

youmengtian

@李东岳感兴趣呀😄

程迪

@李东岳

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

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

yhdthu

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

李东岳

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

程迪

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

程迪

@yhdthu aerodynamics

Him

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

Him

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

浪迹天大

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

wh3296

@程迪
您好，有个问题想请教您。看到您在《OpenFOAM的残差定义》一文中提到了残差定义方式：https://chengdi123000.github.io/2018/01/04/OpenFOAM的残差定义/

关于分母normFactor，有些疑问，不知道为什么要这样定义，有什么数学上的原理吗？
谢谢！