并行效率疑问

程迪

ilu0没有填入，of玩的ldu全是用的不填入的算法，不然ldu数据结构包不住，得像aeroFoam那样玩更高级的数据结构，
ldu矩阵renumber不降低单次迭代的flops数量!那么省出来的时间必须有个说法，要么谱半径更小，收敛所需迭代数减少，要么平均每个flops用时减少(cache命中率提高)。总时间=迭代数每次迭代的flops数每个flops平均用时，对吧。
renumber有效性说明现代cfd程序更多地受限于储存速度，cache命中率啥的，这个东西是有瓶颈的，未来cfd程序的套路可能就得变了。
你说的半带宽那套是针对banded structure的老的数据结构了，算法和数据结构是配套的，单看算法不看数据结构没有任何意义。

李东岳

ilu0没有填入，of玩的ldu全是用的不填入的算法
你说的半带宽那套是针对banded structure的老的数据结构了

在另一个帖子我也看到你说的LDU，这只是OpenFOAM存储矩阵的方式，没有任何算法的存在意义，OpenFOAM编程可以酸则LDU，也可以选择Coordinate list (COO)，也可以选择List of lists (LIL)。ILU0是ILU0算法，LDU是OpenFOAM的存储矩阵的方式。我觉得你弄混了。

另外，矩阵存在带宽是客观存在的，不管你怎么存储。是LDU还是LIL。

矩阵存储和带宽是两个概念：https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_.28CSR.2C_CRS_or_Yale_format.29

对ILU矩阵进行reordering是另一个概念：
Saad, Yousef. Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics, 2003.

程迪

@李东岳

对ILU矩阵进行reordering是另一个概念：

你弄混了，从来没有ILU矩阵，是ldu矩阵吧。

COO/LIL都是sparse matrix的储存格式之一，但是不同储存格式做matrix product或者其他矩阵操作时效率是不一样的，比如CSR和CSC一个做M*x效率高，另一个做transpose(M)*x效率更高。这就是我说的具体数据结构和算法是相关的，一个操作在COO和LDU上的时间复杂度是不一样的。但是很容看出来，ldu矩阵下，很多操作的浮点操作数flops（差不多就是时间复杂度）不随reordering变化。

同时OF默认实现的只有LDU矩阵，没有coo/lil格式的，除了一些GPU扩展。所以OF没法直接调用petsc，也没有依赖别的线性代数软件包，blas，lapack之类的一个没用上，代码也很臃肿。

ILU0是fill in=0的ILU算法，可以基于任何矩阵实现，但是对于sparse matrix来说，fill in为0意味着前后indexing的部分可以不用动，只要把data部分改了就行，甚至指向另一个data_ilu就可以了。OF LDU matrix相当于所有OF里的矩阵的indexing部分都是一样，你能改变的只有data_ptr，所以OF只能玩ilu0，玩ilu1都不行。虽然OF的ldu是很紧凑高效的矩阵内存结构，但是还是有局限性的。

我关注的问题是为什么矩阵减小带宽会带来加速？用半带宽解释我觉得潜在假设就是：同样的矩阵内容，半带宽小==矩阵操作快。但是操作数不变的情况下操作加快了，就是cache的问题了。

李东岳

我关注的问题是为什么矩阵减小带宽会带来加速？

如果去看renumberMesh或者ANSYS关于reorder的描述，并没有说加快计算，只是说减少带宽和提高内存使用率。

Renumbers the cell list in order to reduce the bandwidth, reading and renumbering all fields from all the time directories.

Reordering the domain can improve the computational performance of the solver by rearranging the nodes, faces, and cells in memory. The Mesh/Reorder submenu contains commands for reordering the domain and zones, and also for printing the bandwidth of the present mesh partitions. The domain can be reordered to increase memory access efficiency, and the zones can be reordered for user-interface convenience. The bandwidth provides insight into the distribution of the cells in the zones and in memory.

程迪

这个说得比较清楚。内存效率提升，但是也没说计算浮点操作数减少。

%pylab

Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

A矩阵是对角占优的矩阵

A=array([5 ,-2 ,0 ,0 ,0,
        -2 ,5,-2,0,0,
        0 , -2, 5,-2,0,
        0,0,-2,5,-2,
        0,0,0,-2,5])
A = A.reshape([5,5])
A

array([[ 5, -2,  0,  0,  0],
       [-2,  5, -2,  0,  0],
       [ 0, -2,  5, -2,  0],
       [ 0,  0, -2,  5, -2],
       [ 0,  0,  0, -2,  5]])

P矩阵是排列矩阵，交换第一个元素和第5个元素

P=array([0 ,0 ,0 ,0 ,1,
        0,1,0,0,0,
        0 , 0, 1,0,0,
        0,0,0,1,0,
        1,0,0,0,0])
P = P.reshape([5,5])
P

array([[0, 0, 0, 0, 1],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [1, 0, 0, 0, 0]])

排列矩阵的转置是自身的逆

P.T.dot(P)

array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1]])

作用到矩阵A上，生成矩阵B

B=P.T.dot(A).dot(P)
B

array([[ 5,  0,  0, -2,  0],
       [ 0,  5, -2,  0, -2],
       [ 0, -2,  5, -2,  0],
       [-2,  0, -2,  5,  0],
       [ 0, -2,  0,  0,  5]])

上下三角矩阵和对角矩阵

#GS iteration
L=tril(B,-1)
L

array([[ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0, -2,  0,  0,  0],
       [-2,  0, -2,  0,  0],
       [ 0, -2,  0,  0,  0]])

D=diag(diag(B))
D

array([[5, 0, 0, 0, 0],
       [0, 5, 0, 0, 0],
       [0, 0, 5, 0, 0],
       [0, 0, 0, 5, 0],
       [0, 0, 0, 0, 5]])

U=triu(B,1)
U

array([[ 0,  0,  0, -2,  0],
       [ 0,  0, -2,  0, -2],
       [ 0,  0,  0, -2,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0]])

M_GS is Gauss-Seidel iter. matrix

M_GS=inv(L+D).dot(U)
M_GS

array([[ 0.   ,  0.   ,  0.   , -0.4  ,  0.   ],
       [ 0.   ,  0.   , -0.4  ,  0.   , -0.4  ],
       [ 0.   ,  0.   , -0.16 , -0.4  , -0.16 ],
       [ 0.   ,  0.   , -0.064, -0.32 , -0.064],
       [ 0.   ,  0.   , -0.16 ,  0.   , -0.16 ]])

spectral radius of GS iter.

spectral_radius_GS=max(abs(linalg.eig(M_GS)[0]))
spectral_radius_GS

0.48000000000000026

M_J is Jacobi iter. matrix

M_J=inv(D).dot(U+L)
M_J

array([[ 0. ,  0. ,  0. , -0.4,  0. ],
       [ 0. ,  0. , -0.4,  0. , -0.4],
       [ 0. , -0.4,  0. , -0.4,  0. ],
       [-0.4,  0. , -0.4,  0. ,  0. ],
       [ 0. , -0.4,  0. ,  0. ,  0. ]])

spectral radius of iter. matrix

spectral_radius_J=max(abs(linalg.eig(M_J)[0]))
spectral_radius_J

0.69282032302755092

Compare A with B

L,D,U=tril(A,-1),diag(diag(A)),triu(A,1)
M_GS=inv(L+D).dot(U)
spectral_radius_GS=max(abs(linalg.eig(M_GS)[0]))
spectral_radius_GS

0.47999999999999987

对比可以看出来，reorder之后迭代矩阵谱半径几乎没变化。说明加速主要来源于内存访问效率的提高。

某种意义上来说，说明CPU没吃饱，上菜速度太慢了。以后pde模拟程序的优化方向大概应该是加快上菜速度。。。

Foamer24

@CFDngu 您好！我也遇到了因renumberMesh而不能计算的问题，推测问题和网格有关。没有边界层时，使用renumberMesh没问题且速度会明显变快；在加密边界层至0.00001m后，使用renumberMesh就会出错，当然不使用renumberMesh，也只能单核计算，并行还是会出错。不知道您最近有没有进展，我的问题和算例文件详细描述在 http://www.cfd-china.com/topic/2548/openfoam中cyclicami周期边界在有边界层时出现问题。

如果可以，我还想咨询下，renumberMesh是否有设置精度的选项？或者openfoam重某个文件可以设置精度？个人推测网格尺寸变小，会在关联cyclic或cyclicAMI的周期边界时，出现截断误差导致误判？

Zephyr

目前看下来，GAMG在网格太多了以后效率会急剧下降，而PCG方法网格数多了以后根本收敛不了，在OF自己的框架内应该是无解的。之前看到有国内国外的研究小组都成功把PETSC编译进OF了，楼主如果有一定的编程功底可以自行研究一下。

Calf.Z-DNS

PCG在几百核并行的时候速度更快，但是根本收敛不了吗？GAMG在这么多核的时候效率又很慢，貌似没有最优解？

Calf.Z-DNS

@random_ran 我目前想联系国家这边的超算，请问你有了解国内的超算速度吗？哪里的会比较推荐？谢谢

Calf.Z-DNS

目前了解到PCG在多核并行计算的时候效率会比GAMG更高，那这个多核有大概的一个参考数值吗？之前看到过是五百核。
另外这二者对于计算的准确度上面有没有优劣性，PCG在网格数太多的情况下真的会收敛不了？
最后，多核的时候用scotch会更好吗？之前在论坛上看过对于流向明确的算例，比如圆管，用simple，在流动方向上多分块，效果也还可以。
我目前在400核并行计算7500W的网格，圆管加热流动，用scotch和simple算起来速度差不多。GAMG+DIC跟PCG+DIC的速度也差不多。

不知有何高见，谢谢！

李东岳

@Calf-Z-DNS 我一直对稀疏线性系统感兴趣但是好像有生之年不可能用几年时间去研究他们了。哎。PBiCGStab和PBiCG的算法就在那，但是也没透彻分析过。只不过从个人实践来看我只用PBiCGStab没用过PBiCG。精度没对比过。个人觉得线性系统的精度不足以抹杀掉CFD模型和迭代的精度。你若有结果，欢迎反馈

你后来那个做怎么样了？DNS浮力驱动刘的

yhdthu

fumiya的博客：

https://caefn.com/openfoam/utilities-renumbermesh

yhdthu

@random_ran 请问找到进一步提高效率的方法了么？

李东岳

@yhdthu 老铁最近玩什么呢

yhdthu

@东岳前辈好，最近在算les空化模拟，也卡在计算效率太低上了TAT

yhdthu

所以最终结论是PCG+DIC效率最高是么？
PCG+GAMG是不是稳定性好一些？

刘雄国

OpenFOAM 矩阵计算效率是可以更高,这点没有疑问. 矩阵计算效率高不高主要取决于 MPI 通信
1000 processors 以上不同节点上计算 OpenFOAM 要重构

效率高不高除了算法,计算问题类型, 代码优化等有关系

大牛说下一代cfd 需要大规模计算机网格 + Algorithm + Coding 整体考虑我觉得很有道理

问过openfoam developer, 他的意思是 OpenFOAM is a general CFD solver. They want to solve specific problem to make money (ESI).

刘雄国

@刘雄国在并行效率疑问中说：

大牛说下一代cfd 需要大规模计算机网格 + Algorithm + Coding 整体考虑我觉得很有道理

我说错了, 下一代cfd 需要大规模计算机网格 + CFD Algorithm + CFD Coding + physical网格整体考虑

Micro

@random_ran 在并行效率疑问中说：

以后多核计算前，我都会renumberMesh 的

您好！请问，在进行多核计算时，我应该先decompose,再renumber,还是应该先renumber再decompose呢？谢谢！

浪迹天大

@Micro 网上见到的例子都是先分块，再并行地renumberMesh

CFD中文网

并行效率疑问