时间：2022-04-12 来源：转载于：技术邻——神工坊
在工业仿真领域，对各种现实世界的问题进行数值模拟时，如流体动力学分析、电磁场仿真、结构力学应力应变分析等，其控制方程通常是偏微分方程组，在经过不同方法的隐式离散之后最终都可转化为大型稀疏线性方程组。随着人们对计算精度要求的不断提高，方程组的阶数也从上千阶、几十万阶提高到百万、千万阶甚至更高，所需的计算量以及存储需求也随之迅速膨胀。根据一般经验，方程组求解时间会占总计算时间的70%以上，往往是整个计算过程中的性能瓶颈。如果说求解器是工业CAE软件的核心模块，那么大型稀疏线性方程组的求解技术将毫无疑问是底层求解器的核心。
NASA翼型网格经过离散得到的稀疏矩阵（素材来源于网络）
对于不同类型稀疏矩阵表现差异较大，存在收敛性与收敛速度问题，催生了许多预处理技术（Preconditioners）； 对原矩阵的编辑很少，SpMV(Sparse matrix-vector multiplication)是其核心运算； 内存需求小，求解速度较快，算法复杂度低； 较易实现并行化。
通用、稳定；通过前后处理，能够保证计算的收敛性与精度； 对原矩阵的编辑多（分解、排序、缩放等）； 内存需求大，求解速度慢，算法复杂度更高； 并行度有限。
2014-2021年HPCG性能评测结果对比（素材来源于网络）
直接法的基础是矩阵的分解，常见的分解形式有LU分解、Cholesky分解、QR分解等。稀疏线性方程组的两类常见直接求解算法分别为超节点（Supernodal）方法和多波前（Multifrontal）法，其主要思想是将完整的稀疏矩阵的分解任务转化成许多个相对稠密的子矩阵的分解任务，任务间的依赖关系由消去树（Elimination tree）或其他类似的数据结构来确定。直接法的求解步骤通常分为矩阵重排、符号分解、数值分解与回代求解四个部分。
一个稀疏矩阵与其对应的消去树（来自文献6）
当前，国产超级计算机的峰值性能已达每秒十亿亿次量级，不久便将进入百亿亿次（E级）计算时代，我国的神威E级计算机和天河E级计算机已经蓄势待发。这些国际领先的超级计算机为我国科学与工程计算应用迈进超大规模计算时代、实现更高精细度的数值模拟提供了强力支撑。然而，超大规模计算也给高实用性与高性能的大型稀疏线性方程组求解的算法设计与优化带来了巨大挑战。
1 高 实 用
以SiP封装芯片的电磁-热-力耦合数值模拟为例，其稀疏矩阵具有明显的病态特征（来自文献7）
2 高 性 能
解决问题时间与超级计算机性能趋势对比
对于大规模稀疏线性方程组，原有串行和小规模并行模式下的数据结构和算法容易导致并行求解性能低下或失败。在分布式并行层面，需要解决以下几个问题：一是在数据和任务分解方面，如何设计良好的负载均衡策略、稀疏矩阵的高效存储格式以及计算通信重叠等优化策略；二是在负载均衡的前提下，如何设计以尽力避免节点间的通信；三是在内在串行特性导致并行化困难的算法方面，如何改进数据的分布方式以增加并行性。
左图：AMD霄龙CPU，64核128线程， 右图：英伟达Hopper架构GPU，1.8万核心
我们的探索——UNAP
为了解决高实用性的挑战，UNAP结合异构众核处理器多级并行的特点和稀疏矩阵迭代解法的需求，初步探索了各种预处理方法在众核异构平台上的并行实现技术。UNAP已实现PCG、PBiCGStab、GMRES等Krylov子空间方法以及AMG代数多重网格求解器，预处理器包含Jacobi、DIC、DILU等，未来还将开发直接求解模块，以满足来自各领域的复杂应用需求。
为了达到高性能的目标，UNAP根据国产超算的异构特点，结合其处理器的多级缓冲区，实现了计算/通信混合的并行迭代算法；由于迭代算法的并发度天然较高，UNAP主要通过算法调整、浮点运算替换内存访问、通信同步等手段[9]，充分利用了申威主从众核浮点计算性能高而带宽受限的特点，且通过增加计算比例降低了全局集合通信代价。在共享内存层面，UNAP还可以调用太湖之光超级计算机上的加速工具套件UNAT和向量计算加速库swArrays，充分发挥出从核阵列的计算能力，达到进一步的性能提升。
发动机燃烧室大涡模拟：在某航空发动机全环燃烧室的大涡模拟中，网格量达10亿，采用UNAP作为核心求解模块后，最终的并行规模达到1.6万进程，稀疏线性方程组求解部分加速达到20倍，收敛速度显著提升。 船舶水动力学应用：在某水动力学软件中，通过UNAP的接入（主要使用了求解压力和压力修正方程的 代数多重网格算法和求解速度等方程的预条件稳定双共轭梯度法），获得了在代数求解过程中的自动多级并行能力，在神威·太湖之光超级计算机上完成了千万级网格对标算例计算，并行规模达到万核级别，相对百进程并行效率不低于50%，经测试与商业软件FLUENT相当。 电机设备电磁场分析：在某公司核心求解器向神威平台的移植部署中，采用UNAP代替原有直接求解库，进行了网格量为千万级的电机模型有限元算例并行计算测试，并对比了替换前后算例的节点磁通密度计算结果，UNAP表现良好。
算例 ：方腔驱动流，不可压，上壁面滑移速度为1，其他为固壁边界，方腔的边长为1，雷诺数100。 并行计算结果：以4进程为例，下图给出了方腔在各个进程的分割情况，流线图中可以清楚看到主涡和二次涡。
方腔驱动流的流线图
强扩展性测试：
弱扩展性测试：
[1] Saad Y. Iterative methods for sparse linear systems[M]. Society for Industrial and Applied Mathematics, 2003.
[2] Davis T A, Rajamanickam S, Sid-Lakhdar W M. A survey of direct methods for sparse linear systems[J]. Acta Numerica, 2016, 25: 383-566.
[3] Barrett R, Berry M, Chan T F, et al. Templates for the solution of linear systems: building blocks for iterative methods[M]. Society for Industrial and Applied Mathematics, 1994.
[4] Marjanović V, Gracia J, Glass C W. Performance modeling of the HPCG benchmark[C]//International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems. Springer, Cham, 2014: 172-192.
[5] Cipra B A. The best of the 20th century: Editors name top 10 algorithms[J]. SIAM news, 2000, 33(4): 1-2.
[6] Gupta A, Karypis G, Kumar V. Highly scalable parallel algorithms for sparse matrix factorization[J]. IEEE Transactions on Parallel and Distributed systems, 1997, 8(5): 502-520.
[7] Wang W, Liu Y, Zhao Z, et al. Parallel Multiphysics Simulation of Package Systems Using an Efficient Domain Decomposition Method[J]. Electronics, 2021, 10(2): 158.
[8] 刘伟峰. 高可扩展, 高性能和高实用的稀疏矩阵计算研究进展与挑战[J]. 数值计算与计算机应用, 2020, 41(4): 259.
[9] GU H, REN H U, LIU C, et al. An optimized Chebyshev smoother In GAMG solver of openfoam on sunway Taihulight supercomputer[C]//The 13th OpenFOAM Workshop. 2018.
End
来源：转载于：技术邻——神工坊
