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，不过也分点隐、线隐之类的玩意儿，看上去整得挺高级。
-
可以，但是精度个人感觉取决于对角占优的程度。原方程是：
$$(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实现来看需要对矩阵系数进行操作，里面数据依赖关系决定了可能还是得把矩阵组装出来。
- 如果不用稀疏矩阵，每个标量方程的系数就是
-
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$足够小。我测试了一下，确实相比原矩阵要小。然后我尝试计算了以下矩阵，结果看起来也还
凑合？
上图为通过逆的方法直接求解
上图为LU-SGS非迭代求解咱们先从LU-SGS求解器开始，暂且不讨论效率和线性化，各个击破:sunglasses:
如果我的理解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应该问题不大。 -
@youmengtian
明显Jacobian Free和matrix free不是一回事儿吧。JFNK虽然是Jacobian free，但是Krylov没有preconditioner不是很好使，加preconditioner又很难避免要存、算矩阵，至少算矩阵系数是免不了的，不过这个preconditioner反正是approximation，似乎不需要非常consistent。而且preconditioner整的不好就是并行的性能瓶颈。
-
@youmengtian
Jacobian free是不算Jacobian，JFNK的方法压根就没出现Jacobian，既不存，也不算；Matrix free主要是不存系数矩阵，至于这个矩阵是不是Jacobian无所谓的，但是一般多少都有点儿关系。Matrix free不代表不算矩阵元素，preconditioner免不了要算全部或者部分矩阵元素的，但是往往算完之后可以扔掉一些，不用都存下来。SGS就是一例，按照LU-SGS的搞法，就是Matrix free的，计算过程中用到了某些系数，就算一下，用完了就扔了，下次要用再算。 -
@youmengtian
没，网上有个那个Chun Shen的代码，AETK，我正在看。
https://github.com/chengdi123000/AETKv1
