背景

在工業(yè)仿真領(lǐng)域，對各種現(xiàn)實(shí)世界的問題進(jìn)行數(shù)值模擬時(shí)，如流體動力學(xué)分析、電磁場仿真、結(jié)構(gòu)力學(xué)應(yīng)力應(yīng)變分析等，其控制方程通常是偏微分方程組，在經(jīng)過不同方法的隱式離散之后最終都可轉(zhuǎn)化為大型稀疏線性方程組。隨著人們對計(jì)算精度要求的不斷提高，方程組的階數(shù)也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計(jì)算量以及存儲需求也隨之迅速膨脹。根據(jù)一般經(jīng)驗(yàn)，方程組求解時(shí)間會占總計(jì)算時(shí)間的70%以上，往往是整個(gè)計(jì)算過程中的性能瓶頸。如果說求解器是工業(yè)CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術(shù)將毫無疑問是底層求解器的核心。     

    NASA翼型網(wǎng)格經(jīng)過離散得到的稀疏矩陣（素材來源于網(wǎng)絡(luò)）      

 

 

眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優(yōu)劣，特點(diǎn)比較如下：

迭代法[1]：    

1. 對于不同類型稀疏矩陣表現(xiàn)差異較大，存在收斂性與收斂速度問題，催生了許多預(yù)處理技術(shù)（Preconditioners）；
2. 對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運(yùn)算；
3. 內(nèi)存需求小，求解速度較快，算法復(fù)雜度低；
4. 較易實(shí)現(xiàn)并行化。

直接法[2]：    

1. 通用、穩(wěn)定；通過前后處理，能夠保證計(jì)算的收斂性與精度；
2. 對原矩陣的編輯多（分解、排序、縮放等）；
3. 內(nèi)存需求大，求解速度慢，算法復(fù)雜度更高；
4. 并行度有限。

其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。經(jīng)典的定常迭代法有Jacobi迭代法、Gauss-Seidel迭代法、SOR迭代法等，這些方法均可基于矩陣分裂推導(dǎo)得到；而在數(shù)值模擬中，非定常迭代法則顯得更加重要，常見的有共軛梯度法（Conjugate Gradient CG）、廣義最小殘量方法（Generalized Minimal Residual GMRES）、穩(wěn)定雙共軛梯度法（Biconjugate Gradient Stabilized Bi-CGSTAB）等，這些算法都屬于Krylov子空間方法。 如今高性能共軛梯度（HPCG）測試包已成為國際上評測超級計(jì)算機(jī)性能的主要工具[4]，Krylov子空間方法也被評為20世紀(jì)最偉大的十大算法之一[5]。 還有一些迭代算法則比較特殊，它們基于特殊結(jié)構(gòu)和性質(zhì)，例如我們耳熟能詳?shù)?span style="background-position: initial;background-size: initial;background-repeat: initial;background-attachment: initial;background-origin: initial;background-clip: initial;">多重網(wǎng)格（Multigrid MG）方法，它更多地被稱為數(shù)值計(jì)算領(lǐng)域中一種加速迭代收斂的技術(shù)，而不僅僅是一種單純的算法。    

   

   

2014-2021年HPCG性能評測結(jié)果對比（素材來源于網(wǎng)絡(luò)）     

   

直接法的基礎(chǔ)是矩陣的分解，常見的分解形式有LU分解、Cholesky分解、QR分解等。稀疏線性方程組的兩類常見直接求解算法分別為超節(jié)點(diǎn)（Supernodal）方法和多波前（Multifrontal）法，其主要思想是將完整的稀疏矩陣的分解任務(wù)轉(zhuǎn)化成許多個(gè)相對稠密的子矩陣的分解任務(wù)，任務(wù)間的依賴關(guān)系由消去樹（Elimination tree）或其他類似的數(shù)據(jù)結(jié)構(gòu)來確定。直接法的求解步驟通常分為矩陣重排、符號分解、數(shù)值分解與回代求解四個(gè)部分。    

一個(gè)稀疏矩陣與其對應(yīng)的消去樹（來自文獻(xiàn)6）    

   

 

挑戰(zhàn)

當(dāng)前，國產(chǎn)超級計(jì)算機(jī)的峰值性能已達(dá)每秒十億億次量級，不久便將進(jìn)入百億億次（E級）計(jì)算時(shí)代，我國的神威E級計(jì)算機(jī)和天河E級計(jì)算機(jī)已經(jīng)蓄勢待發(fā)。這些國際領(lǐng)先的超級計(jì)算機(jī)為我國科學(xué)與工程計(jì)算應(yīng)用邁進(jìn)超大規(guī)模計(jì)算時(shí)代、實(shí)現(xiàn)更高精細(xì)度的數(shù)值模擬提供了強(qiáng)力支撐。然而，超大規(guī)模計(jì)算也給高實(shí)用性與高性能的大型稀疏線性方程組求解的算法設(shè)計(jì)與優(yōu)化帶來了巨大挑戰(zhàn)。

對于大規(guī)模稀疏線性方程組，原有串行和小規(guī)模并行模式下的數(shù)據(jù)結(jié)構(gòu)和算法容易導(dǎo)致并行求解性能低下或失敗。在分布式并行層面，需要解決以下幾個(gè)問題：一是在數(shù)據(jù)和任務(wù)分解方面，如何設(shè)計(jì)良好的負(fù)載均衡策略、稀疏矩陣的高效存儲格式以及計(jì)算通信重疊等優(yōu)化策略；二是在負(fù)載均衡的前提下，如何設(shè)計(jì)以盡力避免節(jié)點(diǎn)間的通信；三是在內(nèi)在串行特性導(dǎo)致并行化困難的算法方面，如何改進(jìn)數(shù)據(jù)的分布方式以增加并行性。

在共享內(nèi)存環(huán)境中，稀疏線性方程組求解算法的可擴(kuò)展性問題也需要特別關(guān)注。因?yàn)楝F(xiàn)代多核/眾核處理器上的核數(shù)在可預(yù)見的未來也將越多，在單個(gè)CPU上封裝數(shù)十甚至上百個(gè)有較強(qiáng)處理能力的核心，或是在GPU上封裝成千上萬個(gè)輕量級處理單元將變得非常普遍。如何在這種共享內(nèi)存節(jié)點(diǎn)上實(shí)現(xiàn)細(xì)粒度的并行仍然是很有挑戰(zhàn)的研究內(nèi)容。  

       

    左圖：AMD霄龍CPU，64核128線程， 右圖：英偉達(dá)Hopper架構(gòu)GPU，1.8萬核心   

       

     

     

下面通過幾個(gè)案例說明UNAP的應(yīng)用情況：

• 發(fā)動機(jī)燃燒室大渦模擬： 在某航空發(fā)動機(jī)全環(huán)燃燒室的大渦模擬中，網(wǎng)格量達(dá)10億，采用UNAP作為核心求解模塊后，最終的并行規(guī)模達(dá)到1.6萬進(jìn)程，稀疏線性方程組求解部分加速達(dá)到20倍，收斂速度顯著提升。
• 船舶水動力學(xué)應(yīng)用： 在某水動力學(xué)軟件中，通過UNAP的接入（主要使用了求解壓力和壓力修正方程的   代數(shù)多重網(wǎng)格算法 和求解速度等方程 的預(yù)條件穩(wěn)定雙共軛梯度法 ），獲得了在代數(shù)求解過程中的自動多級并行能力，在神威·太湖之光超級計(jì)算機(jī)上完成了千萬級網(wǎng)格對標(biāo)算例計(jì)算， 并行規(guī)模達(dá)到萬核級別，相對百進(jìn)程 并行效率不低于50% ，經(jīng)測試與商業(yè)軟件FLUENT相當(dāng)。
• 電機(jī)設(shè)備電磁場分析： 在某公司核心求解器向神威平臺的移植部署中，采用UNAP代替原有直接求解庫，進(jìn)行了網(wǎng)格量為千萬級的電機(jī)模型有限元算例并行計(jì)算測試，并對比了替換前后算例的節(jié)點(diǎn)磁通密度計(jì)算結(jié)果，UNAP表現(xiàn)良好。

最后通過一個(gè)典型算例展示UNAP的功能與性能：        

• 算例 ： 方腔驅(qū)動流，不可壓，上壁面滑移速度為1，其他為固壁邊界，方腔的邊長為1，雷諾數(shù)100。
• 并行計(jì)算結(jié)果： 以4進(jìn)程為例，下圖給出了方腔在各個(gè)進(jìn)程的分割情況，流線圖中可以清楚看到主渦和二次渦。

     

方腔驅(qū)動流的流線圖      

   

   

     

下面進(jìn) 行簡單的并行計(jì)算效率測試：                    

• 強(qiáng)擴(kuò)展性測試：

　　固定算例網(wǎng)格整體求解規(guī)模，通過比較不同并行核數(shù)下的計(jì)算時(shí)間可獲得 強(qiáng)擴(kuò)展性 并行計(jì)算效率。本測試采用的網(wǎng)格規(guī)模數(shù)量為2千萬，網(wǎng)格類型為六面體，進(jìn)程數(shù)量分別為64、256、1024，并行核數(shù)規(guī)模分別為4160、16640和66560，統(tǒng)計(jì)時(shí)間為前20步計(jì)算時(shí)間。          

• 弱擴(kuò)展性測試：

　     固定每個(gè)進(jìn)程上的網(wǎng)格規(guī)模，通過比較不同并行核數(shù)下的計(jì)算時(shí)間來獲得 弱擴(kuò)展性 并行計(jì)算效率。本測試采用的網(wǎng)格規(guī)模數(shù)量分別為2千萬、4千萬、8千萬，網(wǎng)格類型為六面體，進(jìn)程數(shù)量分別為64進(jìn)程、128進(jìn)程和256進(jìn)程，并行核數(shù)規(guī)模分別為4160、16640和66560，單個(gè)進(jìn)程網(wǎng)格數(shù)量約為31.25萬，統(tǒng)計(jì)時(shí)間為前20步計(jì)算時(shí)間。（本文作者：趙程鵬）

       

參考文獻(xiàn)

[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] 劉偉峰. 高可擴(kuò)展, 高性能和高實(shí)用的稀疏矩陣計(jì)算研究進(jìn)展與挑戰(zhàn)[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用, 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.

十四五期間，工業(yè)數(shù)字化將是工業(yè)轉(zhuǎn)型升級的主路線?！吧窆し弧北帧八懔x能、協(xié)同創(chuàng)新”的理念，爭做“先進(jìn)算力到仿真算能的轉(zhuǎn)換器”、“離散機(jī)理和垂直仿真場景的連接器”，助力我國工程仿真技術(shù)實(shí)現(xiàn)跨越發(fā)展，支撐重大裝備研制創(chuàng)新和工業(yè)設(shè)計(jì)研發(fā)數(shù)字化轉(zhuǎn)型。