鄒軍,楊帥,程啟問
(清華大學(xué)電機(jī)系, 北京 100084)
構(gòu)建以新能源為主體的新型電力系統(tǒng)是實(shí)現(xiàn)中國(guó)“30·60”“雙碳”戰(zhàn)略目標(biāo)的基礎(chǔ)支撐。新型電力系統(tǒng)包含交流線路、直流線路、高壓電力設(shè)備和大量新能源接入的電力電子裝置。由于新型電力設(shè)備和電力電子裝置的接入,電力系統(tǒng)中將存在大量時(shí)變的和非線性的“元件”,電力系統(tǒng)運(yùn)行特性變得更為復(fù)雜[1-3]。為模擬電力系統(tǒng)的特性,亟須逐步從器件級(jí)、變換器級(jí)、系統(tǒng)級(jí)開展不同尺度建模和仿真,這些仿真的基礎(chǔ)之一是對(duì)高壓電氣設(shè)備和裝置的電磁行為特性進(jìn)行高精度建模。
有限元方法是一種通用的偏微分方程求解技術(shù),其對(duì)于任意幾何形狀和材料都能夠穩(wěn)定獲得數(shù)值解,為實(shí)際工程實(shí)踐提供定量化的指導(dǎo)[4]??茖W(xué)和技術(shù)領(lǐng)域中對(duì)電磁建模的需求不斷增長(zhǎng),催生了以有限元方法為基礎(chǔ)的商業(yè)軟件。商用電磁仿真軟件具有友好的人機(jī)接口,可與計(jì)算機(jī)輔助設(shè)計(jì)(computer aided design,CAD和計(jì)算機(jī)輔助工程computer aided engineering,CAE)的軟件結(jié)合,在實(shí)踐中獲得了巨大的成功。在電磁仿真領(lǐng)域,借助各種變換、因子分解和對(duì)角化技術(shù),求解數(shù)十億自由度的稠密矩陣已經(jīng)成為可能。新型寬帶求解器、等效方法和區(qū)域分解法等計(jì)算技術(shù)可以對(duì)多尺度問題進(jìn)行全波分析,而不必對(duì)實(shí)際問題做過多的簡(jiǎn)化。有人認(rèn)為數(shù)值計(jì)算機(jī)技術(shù)已經(jīng)“完結(jié)”,無需研究新方法和技術(shù),所有的電磁問題可以采用商業(yè)軟件和集成更高算力的計(jì)算機(jī)/工作站最終得以解決。但是,由于新應(yīng)用場(chǎng)景不斷出現(xiàn),例如納米電磁學(xué)、手征材料和強(qiáng)非線性多物理耦合等場(chǎng)合,商業(yè)軟件功能局限性日益突出,對(duì)于計(jì)算規(guī)模在數(shù)千萬單元、瞬態(tài)、非線性和多物理場(chǎng)耦合問題,仿真能力和手段還非常有限。對(duì)于所仿真的系統(tǒng)往往通過增加硬件開銷的方式采用蠻力進(jìn)行計(jì)算,計(jì)算效率無法保證。因此,對(duì)有限元方法和加速技術(shù)的研究重新被學(xué)術(shù)和工業(yè)界人士所重視。
本文聚焦于有限元方法,從高性能有限元技術(shù)的發(fā)展路徑出發(fā),討論有限元技術(shù)發(fā)展的方向和當(dāng)前的瓶頸問題。同時(shí),結(jié)合電力系統(tǒng)電磁仿真的特點(diǎn),簡(jiǎn)要討論了計(jì)算電磁學(xué)的關(guān)鍵技術(shù)問題和可能的發(fā)展路線圖。
在20世紀(jì)初,不同研究領(lǐng)域的科學(xué)家都曾獨(dú)立提出了對(duì)微分方程的變分形式進(jìn)行離散的求解思路[5],其中德裔美國(guó)數(shù)學(xué)家Courant的工作被認(rèn)為是有限元方法的起源:Courant在分析軸的扭轉(zhuǎn)問題時(shí)將軸的橫截面離散成若干三角形單元,并用每個(gè)節(jié)點(diǎn)上未知量的線性插值表示截面上的應(yīng)力分布。在20世紀(jì)50年代初,波音公司的工程師獨(dú)立發(fā)展出了基于能量極小原理的“直接剛性方法”,并廣泛應(yīng)用于機(jī)翼結(jié)構(gòu)的建模仿真[6]。20世紀(jì)年代初,為了突出在有限單元上進(jìn)行積分運(yùn)算的特性,Clough正式提出了“有限元方法”這一名詞。作為Rayleigh-Ritz方法的一種形式,有限元分析開始被學(xué)術(shù)界廣泛接受,不同領(lǐng)域的學(xué)者嘗試在本領(lǐng)域中使用有限元方法,應(yīng)用數(shù)學(xué)家也開始著手建立有限元方法的數(shù)學(xué)理論。1962年至1972年被稱為有限元方法發(fā)展的黃金十年[7]:有限元的求解對(duì)象從橢圓方程拓展到雙曲方程和拋物方程,應(yīng)用場(chǎng)景從結(jié)構(gòu)力學(xué)拓展到熱傳導(dǎo)和流體分析;經(jīng)典的有限元技術(shù)諸如分片測(cè)試、協(xié)調(diào)單元和非協(xié)調(diào)單元、等參變換在這一時(shí)期被提出;針對(duì)有限元方法的誤差分析和收斂性分析理論被建立;在這一時(shí)期也涌現(xiàn)了大量經(jīng)典有限元計(jì)算程序,如NASA發(fā)展的結(jié)構(gòu)有限元分析軟件NASTRAN、John Swanon發(fā)展的通用有限元分析軟件ANSYS等。為了解決傳統(tǒng)變分理論的有限元方法在求解時(shí)出現(xiàn)的數(shù)值振蕩問題,Zienkiewicz等人提出了基于迎風(fēng)格式的Petrov Galerkin有限元方法[8];在此基礎(chǔ)上針對(duì)暫態(tài)問題發(fā)展出了與波傳播特征相結(jié)合的Taylor Galerkin有限元方法[9];1990年前后,Cockburn和舒其望等人提出了具有通量守恒特性的離散Galerkin有限元方法[10]。
隨著計(jì)算電磁學(xué)研究對(duì)象越來越復(fù)雜,模型越來越精細(xì),有限元離散后生成的代數(shù)方程組規(guī)模愈來愈大,這對(duì)現(xiàn)有有限元方法提出了新的挑戰(zhàn)。以高電壓工程中的流注發(fā)展問題為例,受限于CFL條件,迭代的時(shí)間步僅為數(shù)皮秒,實(shí)際仿真中需要進(jìn)行上千步迭代,而每一個(gè)時(shí)間步需要求解一個(gè)涉及上億自由度的問題[11];在仿真計(jì)算高壓絕緣子表面電位時(shí),絕緣子的幾何尺度遠(yuǎn)小于桿塔的幾何尺度,導(dǎo)致整個(gè)計(jì)算區(qū)域離散后產(chǎn)生數(shù)百萬個(gè)節(jié)點(diǎn)[12];大地電阻率反演問題中,采用精細(xì)的像素級(jí)剖分后將產(chǎn)生數(shù)百萬待反演參數(shù),在反演過程中需要使用有限元方法反復(fù)求解正問題[13]。由于特高壓工程的推進(jìn),大量超常規(guī)的高壓設(shè)備也需要重新進(jìn)行設(shè)計(jì):大容量電力電子器件中電-熱-力三種物理場(chǎng)相互耦合,共同決定了開關(guān)性能;針對(duì)高壓套管中電場(chǎng)優(yōu)化需要解決分層介質(zhì)導(dǎo)致的多重尺度問題;在電力變壓器中各類軟磁材料引入的強(qiáng)非線性對(duì)數(shù)值仿真的穩(wěn)定性和精度提出挑戰(zhàn)。總結(jié)來看,電力系統(tǒng)高壓設(shè)備電磁仿真中存在大量三維、時(shí)域以及非線性問題亟待解決,為此,中國(guó)學(xué)者提出了“計(jì)算高壓學(xué)”的概念,而有限元是其核心支撐技術(shù)之一。在新型電力系統(tǒng)中,有大量電力電子設(shè)備的應(yīng)用,考慮到整個(gè)系統(tǒng)的電磁特性,計(jì)算空間尺度可能會(huì)跨越3個(gè)數(shù)量級(jí)以上,目前還不能完全采用場(chǎng)的方法仿真。
高性能計(jì)算,即使用盡可能多的硬件資源完成特定問題的計(jì)算。將一個(gè)串行計(jì)算程序轉(zhuǎn)化成并行程序并移植到高性能計(jì)算機(jī)上運(yùn)行,涉及兩類操作:分解和映射[14]:將計(jì)算任務(wù)分解成一組按照任務(wù)依賴圖并行執(zhí)行的任務(wù),然后映射到各個(gè)計(jì)算單元上執(zhí)行。為了在最短運(yùn)行時(shí)間內(nèi)執(zhí)行完相關(guān)任務(wù),需要根據(jù)實(shí)際計(jì)算單元特點(diǎn)將任務(wù)分解成特定粒度的子任務(wù),同時(shí)最小化單元間通信開銷,減少不同計(jì)算單元之間的交互。由此可見,高性能計(jì)算技術(shù)的應(yīng)用依托于計(jì)算平臺(tái)本身。
在處理器發(fā)展早期,計(jì)算性能的提升主要靠提升主頻和處理器內(nèi)部的隱式并行[15]:提升主頻,即縮短處理器的時(shí)鐘周期,從而提升單位時(shí)間內(nèi)浮點(diǎn)數(shù)運(yùn)算能力;處理器內(nèi)部隱式并行則是通過流水線技術(shù)結(jié)合超長(zhǎng)指令實(shí)現(xiàn)指令級(jí)并行,提高硬件利用率,進(jìn)而提升計(jì)算性能。由于單個(gè)處理器計(jì)算性能低,內(nèi)存小,為了獲得更高算力,需要將大量處理器通過高速通信網(wǎng)絡(luò)連接,構(gòu)成大規(guī)模并行處理機(jī)(mas parallel processor, MPP)[16]。在這一階段,主要是通過區(qū)域分解法將大規(guī)模問題分解成若干小規(guī)模的子問題并分配給各個(gè)結(jié)點(diǎn)進(jìn)行計(jì)算,再通過結(jié)點(diǎn)間的通信傳遞不同子問題間的相互影響,最終實(shí)現(xiàn)對(duì)原問題的求解。區(qū)域分解法本質(zhì)上是通過犧牲計(jì)算時(shí)間,使用有限計(jì)算資源求解大規(guī)模計(jì)算問題,當(dāng)問題規(guī)模較小時(shí),使用區(qū)域分解法將帶來近十倍的耗時(shí)增長(zhǎng)[17]。
隨著芯片制造工藝發(fā)展和指令集的更新,單核處理器性能逐漸逼近理論上限,為進(jìn)一步提升處理器性能,多核(multi-core)構(gòu)架應(yīng)運(yùn)而生:處理器由多個(gè)獨(dú)立的計(jì)算核心構(gòu)成,各個(gè)計(jì)算核心通過共享高速緩存與多線程技術(shù)相結(jié)合,實(shí)現(xiàn)更高計(jì)算性能。多核構(gòu)架受限于高速緩存一致性,隨著核數(shù)增加,緩存信息的維護(hù)成本越來越高,難度越來越大。為進(jìn)一步提升處理器計(jì)算性能,在多核構(gòu)架基礎(chǔ)上發(fā)展出了眾核(many-core)構(gòu)架,即放棄高速緩存一致性,通過堆疊更多計(jì)算核心來提升算力。典型的眾核處理器有通用圖形處理器(GPU)、志強(qiáng)協(xié)處理器(Xeon Phi)等?;诓煌O(shè)計(jì)理念的構(gòu)架具有完全不同的特性,相應(yīng)的處理器適用于處理不同的問題[18]:1)多核處理器:適合順序執(zhí)行復(fù)雜指令,串行計(jì)算性能高;由于線程切換開銷大,并行計(jì)算性能較弱。2)眾核處理器:適合執(zhí)行單指令多數(shù)據(jù)(single instruction multi-data,SIMD)或單指令多線程(single instruction multi-thread,SIMT)操作,并行計(jì)算性能高;由于每個(gè)核上寄存器數(shù)量較少,不適合執(zhí)行復(fù)雜指令。
隨著眾核處理器的出現(xiàn),有限元方法能夠進(jìn)行更細(xì)粒度的分解,實(shí)現(xiàn)單元級(jí)并行,即單個(gè)核心(或單個(gè)線程)處理一個(gè)單元。隨著問題規(guī)模增大,GPU的顯存成為限制求解規(guī)模的主要瓶頸[19]:以2020年最新發(fā)布的A100顯卡為例,其最大顯存僅為40 GB。為解決這一問題,研究人員提出了多GPU聯(lián)合仿真的解決方法。如圖1所示,將整體剛性矩陣分塊,每塊GPU負(fù)責(zé)一塊子矩陣的生成和裝配;求解矩陣方程時(shí),各個(gè)GPU并行完成矩陣-向量乘法,將計(jì)算結(jié)果傳遞給CPU,再由CPU進(jìn)行組裝并分發(fā)給各個(gè)GPU。多GPU聯(lián)合仿真時(shí),由于引入了CPU進(jìn)行收據(jù)收集和分發(fā),需要考慮由此產(chǎn)生的CPU-GPU通信開銷:隨著并行GPU的數(shù)量增加,從CPU向GPU分發(fā)數(shù)據(jù)時(shí)會(huì)造成數(shù)據(jù)總線阻塞,相應(yīng)地通訊開銷會(huì)成倍增長(zhǎng)。如表1—2所示,當(dāng)GPU數(shù)量從2塊增長(zhǎng)為4塊時(shí),從CPU向GPU分發(fā)數(shù)據(jù)的時(shí)間從4 ms增長(zhǎng)為8 ms[20]。
表1 塊GPU聯(lián)合仿真時(shí)各類操作的耗時(shí) (GPU數(shù)量為2)Tab.1 Time-consuming of various operations in block GPU co-simulation (the number of GPU is 2)
圖1 多GPU聯(lián)合仿真模式Fig.1 Multi-GPU co-simulation mode
目前主流服務(wù)器、工作站往往同時(shí)提供多核處理器和眾核處理器,而在上述提到的研究中,計(jì)算任務(wù)完全由GPU執(zhí)行,CPU只起到居中協(xié)調(diào)的作用。為了進(jìn)一步提升算力,需要讓閑置的CPU也參與計(jì)算。為實(shí)現(xiàn)CPU和GPU聯(lián)合的異構(gòu)計(jì)算,需要設(shè)計(jì)合理的任務(wù)分配機(jī)制。研究人員提出了兩類任務(wù)分配策略:靜態(tài)分配——在任務(wù)開始前將任務(wù)分配給各個(gè)計(jì)算單元;動(dòng)態(tài)分配——在任務(wù)執(zhí)行過程中將任務(wù)分配給各個(gè)計(jì)算單元。
表2 塊GPU聯(lián)合仿真時(shí)各類操作耗時(shí) (GPU數(shù)量為4)Tab.2 Time-consuming of various operations in block GPU co-simulation (the number of GPU is 4)
如圖2所示,由于CPU向GPU的通信開銷以及GPU的準(zhǔn)備時(shí)間,在問題規(guī)模較小時(shí)使用CPU執(zhí)行的耗時(shí)要小于GPU執(zhí)行的耗時(shí)?;诖耍墨I(xiàn)[21]制定了靜態(tài)的任務(wù)分配策略:當(dāng)問題規(guī)模小于設(shè)定閾值時(shí),只使用CPU;隨著問題規(guī)模增大,將大部分任務(wù)轉(zhuǎn)移到GPU上執(zhí)行,余下少量任務(wù)留在CPU上執(zhí)行。在特定問題規(guī)模之下,這種靜態(tài)分配策略比僅使用GPU的策略快1.5倍。
圖2 不同問題規(guī)模下CPU和GPU計(jì)算時(shí)間對(duì)比Fig.2 Comparison of CPU and GPU computing time under different problem scales
文獻(xiàn)[22]在利用時(shí)域有限元仿真時(shí)采用了動(dòng)態(tài)任務(wù)分配策略:根據(jù)上一個(gè)時(shí)間步的各個(gè)計(jì)算單元的計(jì)算耗時(shí),調(diào)整下一個(gè)時(shí)間步各個(gè)計(jì)算單元的任務(wù)量,從而實(shí)現(xiàn)動(dòng)態(tài)負(fù)載均衡。最終計(jì)算結(jié)果如表3所示,可以看到在動(dòng)態(tài)任務(wù)分配機(jī)制下,CPU+GPU的異構(gòu)計(jì)算性能優(yōu)于僅使用CPU或GPU的計(jì)算模式。
表3 不同硬件平臺(tái)上,GPU-CPU聯(lián)合仿真加速比Table 3 GPU-CPU co-simulation speedup ratio on different hardware platforms
單個(gè)計(jì)算機(jī)的算力有限,為了進(jìn)一步提升計(jì)算能力,基于多機(jī)互聯(lián)的分布式計(jì)算應(yīng)運(yùn)而生:大量由通信網(wǎng)絡(luò)連接的計(jì)算機(jī)構(gòu)成一張網(wǎng)絡(luò),通過任務(wù)分配調(diào)用網(wǎng)絡(luò)中閑置的計(jì)算資源,實(shí)現(xiàn)巨大算力的使用,這便是網(wǎng)格計(jì)算;而由特殊的高速通信網(wǎng)絡(luò)和同構(gòu)計(jì)算機(jī)互聯(lián)形成的則是計(jì)算集群如常見的超級(jí)計(jì)算機(jī)。隨著虛擬化技術(shù)的發(fā)展,網(wǎng)格計(jì)算進(jìn)一步發(fā)展為如今的云計(jì)算,“計(jì)算”作為一種資源以越來越便捷的方式提供給客戶。
對(duì)于非線性問題和時(shí)域問題,通過線性化、時(shí)域離散后都可化簡(jiǎn)為若干穩(wěn)態(tài)問題的求解。使用有限元方法求解穩(wěn)態(tài)問題可分為3步:?jiǎn)卧仃嚿?、整體矩陣裝配以及線性方程組求解。在實(shí)際計(jì)算中這3步都可能成為瓶頸。如圖3所示,對(duì)于非線性色散介質(zhì)中電場(chǎng)仿真,隨著計(jì)算規(guī)模增長(zhǎng),單元矩陣分析和裝配耗時(shí)占總計(jì)算時(shí)長(zhǎng)的99%[23];而在電磁輻射場(chǎng)仿真中,表4結(jié)果指出線性方程組的求解耗時(shí)在總計(jì)算耗時(shí)中占主導(dǎo)地位[20]。下面,將分別介紹單元矩陣生成、裝配和線性方程組求解這3步的加速技術(shù)。
圖3 不同計(jì)算規(guī)模下整體剛性矩陣生成耗時(shí)在總計(jì)算時(shí)間中占比Fig.3 Proportion of the overall rigid matrix generation time in different calculation scales in the total calculation time
表4 矩陣裝配和求解環(huán)節(jié)的計(jì)算時(shí)間對(duì)比Tab.4 Comparison of calculation time between matrix assembly and solution
2.2.1 單元矩陣生成
單元分析時(shí)僅用到各個(gè)單元自身信息,可以進(jìn)行逐單元并行:對(duì)于低階單元,往往使用單個(gè)線程(或單個(gè)核心)來生成對(duì)應(yīng)的單元?jiǎng)傂跃仃嚕?4];對(duì)于高階單元,單元?jiǎng)傂跃仃囈?guī)模較大,需要分配更多計(jì)算資源來生成一個(gè)單元矩陣[25]。由于不同單元在分析過程中都需要讀入相關(guān)單元的節(jié)點(diǎn)信息,大量的內(nèi)存讀寫操作使得這類傳統(tǒng)的并行方法效率提升有限。另一種加速思路是將單元分析過程進(jìn)行張量化改寫。
2.2.2 整體矩陣裝配
獲得各個(gè)單元的局部剛性矩陣后,需要裝配形成整體剛性矩陣。由于一個(gè)節(jié)點(diǎn)的所有鄰接單元矩陣對(duì)這一結(jié)點(diǎn)均有貢獻(xiàn),在裝配時(shí)的歸并操作將導(dǎo)致內(nèi)存的沖突讀寫。為解決這一問題,一種常見的解決思路是染色法:對(duì)計(jì)算區(qū)域內(nèi)單元進(jìn)行染色,使得相鄰單元不同色。每次將具有相同顏色單元的局部剛性矩陣裝配至整體剛性矩陣中,從而避免了競(jìng)爭(zhēng)沖突;另一種解決方法是將矩陣裝配環(huán)節(jié)和后續(xù)的線性方程組求解環(huán)節(jié)相結(jié)合,如逐單元(element by element)有限元[26]、節(jié)點(diǎn)分解(nodal domain decomposition)有限元[27]等,這些免裝配方法實(shí)際上都可視為多波前方法的變種[28]。
2.2.3 稀疏線性方程組求解
通過單元矩陣生成和裝配后,原偏微分方程被離散成大型稀疏線性方程組。如何求解稀疏線性代數(shù)方程組一直是高性能計(jì)算領(lǐng)域的核心議題。針對(duì)大型稀疏線性方程組的解法可以分為兩類。
第一類是基于高斯消去法的直接求解器。直接計(jì)算方法可分為填充元優(yōu)化、符號(hào)分解、數(shù)值分解、消元回代4步。其中填充元優(yōu)化是對(duì)矩陣進(jìn)行重排序,在保證數(shù)值穩(wěn)定性的同時(shí),使得后續(xù)矩陣分解時(shí)產(chǎn)生盡可能少的填充元,從而提高計(jì)算效率。常見的方法如最小度算法、重復(fù)剖分算法以及多重剖分算法。而符號(hào)分解則是現(xiàn)代快速直接解法特有的技術(shù),通過使用消去樹這一數(shù)據(jù)結(jié)構(gòu),確定分解因子矩陣的非零元位置,指導(dǎo)分解路徑,如基于一般消去樹的多波前方法和基于超節(jié)點(diǎn)消去樹的多元最小度方法等。直接解法的優(yōu)點(diǎn)是計(jì)算穩(wěn)定可靠,因此也成為大量通用有限元計(jì)算軟件中的默認(rèn)求解器;其缺點(diǎn)是計(jì)算過程中產(chǎn)生的分解因子、波前矩陣等中間變量規(guī)模不定,隨著問題規(guī)模增大,超出內(nèi)存容量后超出部分只能存儲(chǔ)在硬盤中,硬盤的I/O操作成為計(jì)算瓶頸制約了直接解法的性能。
第二類是迭代解法。預(yù)計(jì)于Gauss消去法的直接解法不同,迭代解法是在特定的解空間中不斷尋找更接近真解的元素,其中最常見的是基于krylov子空間的迭代方法:如先驗(yàn)共軛梯度迭代算法(preconditioned conjugate gradient,PCG)方法、雙共軛梯度(biconjugate gradient,BiCG)方法等。在迭代求解過程中需要不斷進(jìn)行稀疏矩陣向量乘法操作(sparse matrix-vector multiplication,SpMV)來獲得當(dāng)前解對(duì)應(yīng)的殘差,并基于此進(jìn)行修正。如何高效地執(zhí)行SpMV操作成為該類解法的研究重點(diǎn)。研究人員通過調(diào)整稀疏矩陣的存儲(chǔ)格式來實(shí)現(xiàn)內(nèi)存的合并讀寫,優(yōu)化SpMV操作的性能。常見的存儲(chǔ)格式可分為兩類:行主元的COO格式、CSR格式、CTO格式以及列主元的ELL格式,Sliced ELL格式等[29]。針對(duì)硬件配置選擇特定的存儲(chǔ)格式,在存儲(chǔ)空間最小化同時(shí)通過緩存的合并讀寫實(shí)現(xiàn)時(shí)空效率的最大化。
迭代求解器中另一個(gè)關(guān)鍵技術(shù)是預(yù)條件處理。隨著問題規(guī)模增大,剛性矩陣的條件數(shù)越來越大,通過預(yù)條件處理能夠顯著提升迭代算法的收斂性能。在并行求解器中,剛性矩陣往往是分塊存儲(chǔ),如何根據(jù)局部矩陣信息生成預(yù)條件矩陣是并行求解器研究的熱點(diǎn),目前發(fā)展的并行預(yù)條件矩陣有:不完全LU分解預(yù)條件子,分塊Jacobi預(yù)條件子,最小二乘多項(xiàng)式預(yù)條件子以及多重網(wǎng)格預(yù)條件子等。相較于直接求解器,迭代求解器計(jì)算過程中對(duì)存儲(chǔ)的需求固定,借助區(qū)域分解方法容易實(shí)現(xiàn)稀疏矩陣-向量乘法操作的并行化,適合求解較大規(guī)模問題;而迭代解法的缺點(diǎn)也非常明顯:無法保證收斂、生成性態(tài)優(yōu)良的預(yù)條件矩陣較為耗時(shí)。
縱觀有限元發(fā)展歷史,求解器經(jīng)歷了迭代方法(Gauss-Seidel方法)、直接方法(變帶寬方法)、PCG方法再到多波前等直接解法的螺旋式發(fā)展[30]。如今直接解法能夠有效求解百萬自由度的有限元方程,而迭代解法多用于更大規(guī)模問題[31]。無論是直接解法還是迭代解法,性能的提升都是通過發(fā)掘稀疏矩陣自身結(jié)構(gòu)并與計(jì)算機(jī)硬件相結(jié)合來實(shí)現(xiàn),隨著計(jì)算機(jī)硬件發(fā)展,如何利用多核和眾核處理器實(shí)現(xiàn)高效求解成為新的問題。
隨著物聯(lián)網(wǎng)(internet of things,IoT)技術(shù)的成熟,實(shí)時(shí)獲取設(shè)備運(yùn)行過程中的各種參數(shù)成為可能。為了有效利用這些參數(shù)對(duì)設(shè)備全生命周期進(jìn)行預(yù)測(cè)和控制,數(shù)字孿生(digital twin)這一概念被提出。由于實(shí)際設(shè)備的物理模型過于復(fù)雜,傳統(tǒng)數(shù)值仿真往往采用簡(jiǎn)化模型:如使用集總參數(shù)模型、將非線性時(shí)變的材料特性線性化等。仿真模型和實(shí)際模型間存在較大的差異,極大地限制了仿真結(jié)果的準(zhǔn)確性,因此傳統(tǒng)數(shù)值仿真僅作為優(yōu)化設(shè)計(jì)的輔助工具。大量設(shè)備內(nèi)部參數(shù)的出現(xiàn)使得數(shù)值仿真模型能夠深入設(shè)備內(nèi)部,盡可能真實(shí)地還原實(shí)際運(yùn)行工況,在此基礎(chǔ)上得到的仿真計(jì)算結(jié)果準(zhǔn)確度大大提高,基于此能夠進(jìn)行對(duì)設(shè)備運(yùn)行過程中出現(xiàn)的故障進(jìn)行檢測(cè)、定位以及預(yù)測(cè)。
數(shù)字孿生概念的出現(xiàn)極大地拓寬了電磁場(chǎng)數(shù)值仿真的應(yīng)用場(chǎng)景,同時(shí)也提出新的挑戰(zhàn):精細(xì)的模型往往兼具時(shí)域、三維、非線性等特征,某些復(fù)雜的問題還涉及熱場(chǎng)、流場(chǎng)、應(yīng)力場(chǎng)等多物理場(chǎng)耦合。如何利用邊緣計(jì)算、云計(jì)算等新興技術(shù)實(shí)現(xiàn)穩(wěn)定高效地求解是一個(gè)亟待解決的問題。
在國(guó)際科學(xué)計(jì)算社區(qū)中不同研究團(tuán)隊(duì)分別提出并發(fā)布了用于偏微分方程求解的高性能數(shù)值計(jì)算庫:Gmsh提供二維、三維結(jié)構(gòu)的網(wǎng)格剖分;METIS則能夠?qū)W(wǎng)格編號(hào)進(jìn)行重排、分解;針對(duì)大型稀疏代數(shù)方程組,MUMPS提供了基于多波前方法的直接求解器,Hypre、AGMG等軟件則提供了基于代數(shù)多重網(wǎng)格方法的迭代求解器;可移植可拓展科學(xué)計(jì)算工具箱PETSc專門用于求解偏微分方程;Trillinos使用面向?qū)ο蟮能浖蚣?,通過調(diào)用各個(gè)組件實(shí)現(xiàn)大規(guī)模復(fù)雜多物理場(chǎng)求解;另外還存在大量基于上述計(jì)算組件二次開發(fā)的有限元計(jì)算軟件,如Deal.Ⅱ、FENICs、MFEM等。這些開源的高性能計(jì)算軟件是當(dāng)前偏微分方程數(shù)值解研究的重要基礎(chǔ),在計(jì)算流體力學(xué)、計(jì)算力學(xué)等領(lǐng)域都得到廣泛的應(yīng)用,但在計(jì)算電磁學(xué)領(lǐng)域未見報(bào)道。
近年來,大數(shù)據(jù)和人工智能技術(shù)突飛猛進(jìn),深度神經(jīng)網(wǎng)絡(luò)技術(shù)給傳統(tǒng)科學(xué)計(jì)算帶來新的詮釋。神經(jīng)計(jì)算框架之下,科學(xué)計(jì)算問題可分為3類[32]:第1類問題是完全由數(shù)據(jù)驅(qū)動(dòng)的問題,常見的如對(duì)風(fēng)、光等自然資源的實(shí)時(shí)預(yù)測(cè)。由于實(shí)際模型未知,往往直接簡(jiǎn)化為黑箱模型進(jìn)行處理,采用神經(jīng)網(wǎng)絡(luò)對(duì)數(shù)據(jù)進(jìn)行擬合、挖掘;第2類問題則是基于部分物理機(jī)理和少量測(cè)量數(shù)據(jù)的不確定性模型,常見的如逆散射問題,即控制方程已知,根據(jù)測(cè)量數(shù)據(jù)反演出待求的模型參數(shù)。這一類問題往往可以寫為PDE約束下的優(yōu)化問題,借助神經(jīng)網(wǎng)絡(luò)這一優(yōu)化求解器進(jìn)行求解;第3類問題是基于已知模型的確定性問題,經(jīng)典計(jì)算電磁學(xué)問題往往屬于這一類。有限元過程可以等價(jià)改寫為神經(jīng)網(wǎng)絡(luò)的形式(FEA-net),即每個(gè)神經(jīng)元的輸入為單元幾何信息,輸出為單元?jiǎng)傂跃仃嚕?3]。應(yīng)用神經(jīng)網(wǎng)絡(luò)技術(shù)求解上述3類問題具有兩個(gè)顯著優(yōu)勢(shì):一是借助各類演化算法能夠處理復(fù)雜的非線性優(yōu)化問題;二是神經(jīng)網(wǎng)絡(luò)本身是對(duì)計(jì)算任務(wù)的細(xì)粒度分解,能夠?qū)崿F(xiàn)高性能計(jì)算。將經(jīng)典的計(jì)算電磁學(xué)問題與如今蓬勃發(fā)展的神經(jīng)計(jì)算相結(jié)合,也是未來的發(fā)展方向之一。
高性能計(jì)算技術(shù)的發(fā)展與高性能計(jì)算機(jī)的發(fā)展是密切相關(guān)的,基于高性能計(jì)算技術(shù)的有限元方法在計(jì)算電磁學(xué)中的應(yīng)用經(jīng)歷了從最開始只有單核CPU,發(fā)展到多核處理器,再到多核處理器和眾核處理器并存的異構(gòu)計(jì)算環(huán)境。隨著硬件平臺(tái)在不斷更新?lián)Q代,單位計(jì)算資源的成本越來越低:1997年,提供1TFlops計(jì)算能力的硬件需要花費(fèi)4 600萬美元,而在2017年,使用相同的計(jì)算資源僅需要30美元[34]。相較于發(fā)展精巧的計(jì)算方法,高效利用“廉價(jià)”的硬件資源求解復(fù)雜的電磁學(xué)問題更為重要,也更加有效。
在目前高性能有限元計(jì)算方案中,還存在以下問題:其一,現(xiàn)有研究針對(duì)特定硬件的加速技術(shù),這些技術(shù)往往受限于特定的問題或特定的計(jì)算平臺(tái),尚缺乏一款具有可移植性、可拓展性以及較高硬件利用率的計(jì)算程序;其二,除多核與眾核處理器外,張量處理器 TPU[35]、神經(jīng)計(jì)算單元 NPU[36]、量子處理器[37]等新興的計(jì)算單元也在不斷發(fā)展,有限元方法與這些硬件的結(jié)合也是一個(gè)亟待解決的問題。
有限元方法作為一種通用的數(shù)值求解技術(shù),以其強(qiáng)大的計(jì)算能力,在電磁仿真領(lǐng)域獲得巨大成功。由于新應(yīng)用場(chǎng)景的出現(xiàn),電力系統(tǒng)高壓設(shè)備電磁仿真對(duì)有限元方法提出了新的需求。從技術(shù)路徑發(fā)展上,高性能有限元技術(shù)是上述問題的一個(gè)可能的解決方案。本文介紹了高性能有限元方法的主要技術(shù)要點(diǎn),討論其相關(guān)的發(fā)展途徑,希望對(duì)高性能有限元方法在電力系統(tǒng)高壓設(shè)備仿真的應(yīng)用有所推動(dòng)。