鄒軍，楊帥，程啟問(wèn)

（清華大學(xué)電機(jī)系， 北京 100084）

0 引言

構(gòu)建以新能源為主體的新型電力系統(tǒng)是實(shí)現(xiàn)中國(guó)“30·60”“雙碳”戰(zhàn)略目標(biāo)的基礎(chǔ)支撐。新型電力系統(tǒng)包含交流線(xiàn)路、直流線(xiàn)路、高壓電力設(shè)備和大量新能源接入的電力電子裝置。由于新型電力設(shè)備和電力電子裝置的接入，電力系統(tǒng)中將存在大量時(shí)變的和非線(xiàn)性的“元件”，電力系統(tǒng)運(yùn)行特性變得更為復(fù)雜［1-3］。為模擬電力系統(tǒng)的特性，亟須逐步從器件級(jí)、變換器級(jí)、系統(tǒng)級(jí)開(kāi)展不同尺度建模和仿真，這些仿真的基礎(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ì)多尺度問(wèn)題進(jìn)行全波分析，而不必對(duì)實(shí)際問(wèn)題做過(guò)多的簡(jiǎn)化。有人認(rèn)為數(shù)值計(jì)算機(jī)技術(shù)已經(jīng)“完結(jié)”，無(wú)需研究新方法和技術(shù)，所有的電磁問(wèn)題可以采用商業(yè)軟件和集成更高算力的計(jì)算機(jī)/工作站最終得以解決。但是，由于新應(yīng)用場(chǎng)景不斷出現(xiàn)，例如納米電磁學(xué)、手征材料和強(qiáng)非線(xiàn)性多物理耦合等場(chǎng)合，商業(yè)軟件功能局限性日益突出，對(duì)于計(jì)算規(guī)模在數(shù)千萬(wàn)單元、瞬態(tài)、非線(xiàn)性和多物理場(chǎng)耦合問(wèn)題，仿真能力和手段還非常有限。對(duì)于所仿真的系統(tǒng)往往通過(guò)增加硬件開(kāi)銷(xiāo)的方式采用蠻力進(jìn)行計(jì)算，計(jì)算效率無(wú)法保證。因此，對(duì)有限元方法和加速技術(shù)的研究重新被學(xué)術(shù)和工業(yè)界人士所重視。

本文聚焦于有限元方法，從高性能有限元技術(shù)的發(fā)展路徑出發(fā)，討論有限元技術(shù)發(fā)展的方向和當(dāng)前的瓶頸問(wèn)題。同時(shí)，結(jié)合電力系統(tǒng)電磁仿真的特點(diǎn)，簡(jiǎn)要討論了計(jì)算電磁學(xué)的關(guān)鍵技術(shù)問(wèn)題和可能的發(fā)展路線(xiàn)圖。

1 工程電磁場(chǎng)有限元方法回顧與挑戰(zhàn)

在20世紀(jì)初，不同研究領(lǐng)域的科學(xué)家都曾獨(dú)立提出了對(duì)微分方程的變分形式進(jìn)行離散的求解思路［5］，其中德裔美國(guó)數(shù)學(xué)家Courant的工作被認(rèn)為是有限元方法的起源：Courant在分析軸的扭轉(zhuǎn)問(wèn)題時(shí)將軸的橫截面離散成若干三角形單元，并用每個(gè)節(jié)點(diǎn)上未知量的線(xiàn)性插值表示截面上的應(yīng)力分布。在20世紀(jì)50年代初，波音公司的工程師獨(dú)立發(fā)展出了基于能量極小原理的“直接剛性方法”，并廣泛應(yīng)用于機(jī)翼結(jié)構(gòu)的建模仿真［6］。20世紀(jì)年代初，為了突出在有限單元上進(jìn)行積分運(yùn)算的特性，Clough正式提出了“有限元方法”這一名詞。作為Rayleigh-Ritz方法的一種形式，有限元分析開(kāi)始被學(xué)術(shù)界廣泛接受，不同領(lǐng)域的學(xué)者嘗試在本領(lǐng)域中使用有限元方法，應(yīng)用數(shù)學(xué)家也開(kāi)始著手建立有限元方法的數(shù)學(xué)理論。1962年至1972年被稱(chēng)為有限元方法發(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ù)值振蕩問(wèn)題，Zienkiewicz等人提出了基于迎風(fēng)格式的Petrov Galerkin有限元方法［8］；在此基礎(chǔ)上針對(duì)暫態(tài)問(wèn)題發(fā)展出了與波傳播特征相結(jié)合的Taylor Galerkin有限元方法［9］；1990年前后，Cockburn和舒其望等人提出了具有通量守恒特性的離散Galerkin有限元方法［10］。

隨著計(jì)算電磁學(xué)研究對(duì)象越來(lái)越復(fù)雜，模型越來(lái)越精細(xì)，有限元離散后生成的代數(shù)方程組規(guī)模愈來(lái)愈大，這對(duì)現(xiàn)有有限元方法提出了新的挑戰(zhàn)。以高電壓工程中的流注發(fā)展問(wèn)題為例，受限于CFL條件，迭代的時(shí)間步僅為數(shù)皮秒，實(shí)際仿真中需要進(jìn)行上千步迭代，而每一個(gè)時(shí)間步需要求解一個(gè)涉及上億自由度的問(wèn)題［11］；在仿真計(jì)算高壓絕緣子表面電位時(shí)，絕緣子的幾何尺度遠(yuǎn)小于桿塔的幾何尺度，導(dǎo)致整個(gè)計(jì)算區(qū)域離散后產(chǎn)生數(shù)百萬(wàn)個(gè)節(jié)點(diǎn)［12］；大地電阻率反演問(wèn)題中，采用精細(xì)的像素級(jí)剖分后將產(chǎn)生數(shù)百萬(wàn)待反演參數(shù)，在反演過(guò)程中需要使用有限元方法反復(fù)求解正問(wèn)題［13］。由于特高壓工程的推進(jìn)，大量超常規(guī)的高壓設(shè)備也需要重新進(jìn)行設(shè)計(jì)：大容量電力電子器件中電-熱-力三種物理場(chǎng)相互耦合，共同決定了開(kāi)關(guān)性能；針對(duì)高壓套管中電場(chǎng)優(yōu)化需要解決分層介質(zhì)導(dǎo)致的多重尺度問(wèn)題；在電力變壓器中各類(lèi)軟磁材料引入的強(qiáng)非線(xiàn)性對(duì)數(shù)值仿真的穩(wěn)定性和精度提出挑戰(zhàn)?？偨Y(jié)來(lái)看，電力系統(tǒng)高壓設(shè)備電磁仿真中存在大量三維、時(shí)域以及非線(xiàn)性問(wèn)題亟待解決，為此，中國(guó)學(xué)者提出了“計(jì)算高壓學(xué)”的概念，而有限元是其核心支撐技術(shù)之一。在新型電力系統(tǒng)中，有大量電力電子設(shè)備的應(yīng)用，考慮到整個(gè)系統(tǒng)的電磁特性，計(jì)算空間尺度可能會(huì)跨越3個(gè)數(shù)量級(jí)以上，目前還不能完全采用場(chǎng)的方法仿真。

2 高性能計(jì)算有限元方法的應(yīng)用進(jìn)展

2.1 高性能計(jì)算硬件發(fā)展

高性能計(jì)算，即使用盡可能多的硬件資源完成特定問(wèn)題的計(jì)算。將一個(gè)串行計(jì)算程序轉(zhuǎn)化成并行程序并移植到高性能計(jì)算機(jī)上運(yùn)行，涉及兩類(lèi)操作：分解和映射［14］：將計(jì)算任務(wù)分解成一組按照任務(wù)依賴(lài)圖并行執(zhí)行的任務(wù)，然后映射到各個(gè)計(jì)算單元上執(zhí)行。為了在最短運(yùn)行時(shí)間內(nèi)執(zhí)行完相關(guān)任務(wù)，需要根據(jù)實(shí)際計(jì)算單元特點(diǎn)將任務(wù)分解成特定粒度的子任務(wù)，同時(shí)最小化單元間通信開(kāi)銷(xiāo)，減少不同計(jì)算單元之間的交互。由此可見(jiàn)，高性能計(jì)算技術(shù)的應(yīng)用依托于計(jì)算平臺(tái)本身。

在處理器發(fā)展早期，計(jì)算性能的提升主要靠提升主頻和處理器內(nèi)部的隱式并行［15］：提升主頻，即縮短處理器的時(shí)鐘周期，從而提升單位時(shí)間內(nèi)浮點(diǎn)數(shù)運(yùn)算能力；處理器內(nèi)部隱式并行則是通過(guò)流水線(xiàn)技術(shù)結(jié)合超長(zhǎng)指令實(shí)現(xiàn)指令級(jí)并行，提高硬件利用率，進(jìn)而提升計(jì)算性能。由于單個(gè)處理器計(jì)算性能低，內(nèi)存小，為了獲得更高算力，需要將大量處理器通過(guò)高速通信網(wǎng)絡(luò)連接，構(gòu)成大規(guī)模并行處理機(jī)（mas parallel processor， MPP）［16］。在這一階段，主要是通過(guò)區(qū)域分解法將大規(guī)模問(wèn)題分解成若干小規(guī)模的子問(wèn)題并分配給各個(gè)結(jié)點(diǎn)進(jìn)行計(jì)算，再通過(guò)結(jié)點(diǎn)間的通信傳遞不同子問(wèn)題間的相互影響，最終實(shí)現(xiàn)對(duì)原問(wèn)題的求解。區(qū)域分解法本質(zhì)上是通過(guò)犧牲計(jì)算時(shí)間，使用有限計(jì)算資源求解大規(guī)模計(jì)算問(wèn)題，當(dāng)問(wèn)題規(guī)模較小時(shí)，使用區(qū)域分解法將帶來(lái)近十倍的耗時(shí)增長(zhǎng)［17］。

隨著芯片制造工藝發(fā)展和指令集的更新，單核處理器性能逐漸逼近理論上限，為進(jìn)一步提升處理器性能，多核（multi-core）構(gòu)架應(yīng)運(yùn)而生：處理器由多個(gè)獨(dú)立的計(jì)算核心構(gòu)成，各個(gè)計(jì)算核心通過(guò)共享高速緩存與多線(xiàn)程技術(shù)相結(jié)合，實(shí)現(xiàn)更高計(jì)算性能。多核構(gòu)架受限于高速緩存一致性，隨著核數(shù)增加，緩存信息的維護(hù)成本越來(lái)越高，難度越來(lái)越大。為進(jìn)一步提升處理器計(jì)算性能，在多核構(gòu)架基礎(chǔ)上發(fā)展出了眾核（many-core）構(gòu)架，即放棄高速緩存一致性，通過(guò)堆疊更多計(jì)算核心來(lái)提升算力。典型的眾核處理器有通用圖形處理器（GPU）、志強(qiáng)協(xié)處理器（Xeon Phi）等?；诓煌O(shè)計(jì)理念的構(gòu)架具有完全不同的特性，相應(yīng)的處理器適用于處理不同的問(wèn)題［18］：1）多核處理器：適合順序執(zhí)行復(fù)雜指令，串行計(jì)算性能高；由于線(xiàn)程切換開(kāi)銷(xiāo)大，并行計(jì)算性能較弱。2）眾核處理器：適合執(zhí)行單指令多數(shù)據(jù)（single instruction multi-data，SIMD）或單指令多線(xiàn)程（single instruction multi-thread，SIMT）操作，并行計(jì)算性能高；由于每個(gè)核上寄存器數(shù)量較少，不適合執(zhí)行復(fù)雜指令。

隨著眾核處理器的出現(xiàn)，有限元方法能夠進(jìn)行更細(xì)粒度的分解，實(shí)現(xiàn)單元級(jí)并行，即單個(gè)核心（或單個(gè)線(xiàn)程）處理一個(gè)單元。隨著問(wèn)題規(guī)模增大，GPU的顯存成為限制求解規(guī)模的主要瓶頸［19］：以2020年最新發(fā)布的A100顯卡為例，其最大顯存僅為40 GB。為解決這一問(wèn)題，研究人員提出了多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通信開(kāi)銷(xiāo)：隨著并行GPU的數(shù)量增加，從CPU向GPU分發(fā)數(shù)據(jù)時(shí)會(huì)造成數(shù)據(jù)總線(xiàn)阻塞，相應(yīng)地通訊開(kāi)銷(xiāo)會(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í)各類(lèi)操作的耗時(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ī)制。研究人員提出了兩類(lèi)任務(wù)分配策略：靜態(tài)分配——在任務(wù)開(kāi)始前將任務(wù)分配給各個(gè)計(jì)算單元；動(dòng)態(tài)分配——在任務(wù)執(zhí)行過(guò)程中將任務(wù)分配給各個(gè)計(jì)算單元。

表2 塊GPU聯(lián)合仿真時(shí)各類(lèi)操作耗時(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的通信開(kāi)銷(xiāo)以及GPU的準(zhǔn)備時(shí)間，在問(wèn)題規(guī)模較小時(shí)使用CPU執(zhí)行的耗時(shí)要小于GPU執(zhí)行的耗時(shí)?；诖?，文獻(xiàn)［21］制定了靜態(tài)的任務(wù)分配策略：當(dāng)問(wèn)題規(guī)模小于設(shè)定閾值時(shí)，只使用CPU；隨著問(wèn)題規(guī)模增大，將大部分任務(wù)轉(zhuǎn)移到GPU上執(zhí)行，余下少量任務(wù)留在CPU上執(zhí)行。在特定問(wèn)題規(guī)模之下，這種靜態(tài)分配策略比僅使用GPU的策略快1.5倍。

圖2 不同問(wèn)題規(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ò)，通過(guò)任務(wù)分配調(diào)用網(wǎng)絡(luò)中閑置的計(jì)算資源，實(shí)現(xiàn)巨大算力的使用，這便是網(wǎng)格計(jì)算；而由特殊的高速通信網(wǎng)絡(luò)和同構(gòu)計(jì)算機(jī)互聯(lián)形成的則是計(jì)算集群如常見(jiàn)的超級(jí)計(jì)算機(jī)。隨著虛擬化技術(shù)的發(fā)展，網(wǎng)格計(jì)算進(jìn)一步發(fā)展為如今的云計(jì)算，“計(jì)算”作為一種資源以越來(lái)越便捷的方式提供給客戶(hù)。

2.2 有限元方法加速技術(shù)

對(duì)于非線(xiàn)性問(wèn)題和時(shí)域問(wèn)題，通過(guò)線(xiàn)性化、時(shí)域離散后都可化簡(jiǎn)為若干穩(wěn)態(tài)問(wèn)題的求解。使用有限元方法求解穩(wěn)態(tài)問(wèn)題可分為3步：?jiǎn)卧仃嚿?、整體矩陣裝配以及線(xiàn)性方程組求解。在實(shí)際計(jì)算中這3步都可能成為瓶頸。如圖3所示，對(duì)于非線(xiàn)性色散介質(zhì)中電場(chǎng)仿真，隨著計(jì)算規(guī)模增長(zhǎng)，單元矩陣分析和裝配耗時(shí)占總計(jì)算時(shí)長(zhǎng)的99%［23］；而在電磁輻射場(chǎng)仿真中，表4結(jié)果指出線(xiàn)性方程組的求解耗時(shí)在總計(jì)算耗時(shí)中占主導(dǎo)地位［20］。下面，將分別介紹單元矩陣生成、裝配和線(xiàn)性方程組求解這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è)線(xiàn)程（或單個(gè)核心）來(lái)生成對(duì)應(yīng)的單元?jiǎng)傂跃仃嚕?4］；對(duì)于高階單元，單元?jiǎng)傂跃仃囈?guī)模較大，需要分配更多計(jì)算資源來(lái)生成一個(gè)單元矩陣［25］。由于不同單元在分析過(guò)程中都需要讀入相關(guān)單元的節(jié)點(diǎn)信息，大量的內(nèi)存讀寫(xiě)操作使得這類(lèi)傳統(tǒng)的并行方法效率提升有限。另一種加速思路是將單元分析過(guò)程進(jìn)行張量化改寫(xiě)。

2.2.2 整體矩陣裝配

獲得各個(gè)單元的局部剛性矩陣后，需要裝配形成整體剛性矩陣。由于一個(gè)節(jié)點(diǎn)的所有鄰接單元矩陣對(duì)這一結(jié)點(diǎn)均有貢獻(xiàn)，在裝配時(shí)的歸并操作將導(dǎo)致內(nèi)存的沖突讀寫(xiě)。為解決這一問(wèn)題，一種常見(jiàn)的解決思路是染色法：對(duì)計(jì)算區(qū)域內(nèi)單元進(jìn)行染色，使得相鄰單元不同色。每次將具有相同顏色單元的局部剛性矩陣裝配至整體剛性矩陣中，從而避免了競(jìng)爭(zhēng)沖突；另一種解決方法是將矩陣裝配環(huán)節(jié)和后續(xù)的線(xiàn)性方程組求解環(huán)節(jié)相結(jié)合，如逐單元（element by element）有限元［26］、節(jié)點(diǎn)分解（nodal domain decomposition）有限元［27］等，這些免裝配方法實(shí)際上都可視為多波前方法的變種［28］。

2.2.3 稀疏線(xiàn)性方程組求解

通過(guò)單元矩陣生成和裝配后，原偏微分方程被離散成大型稀疏線(xiàn)性方程組。如何求解稀疏線(xiàn)性代數(shù)方程組一直是高性能計(jì)算領(lǐng)域的核心議題。針對(duì)大型稀疏線(xiàn)性方程組的解法可以分為兩類(lèi)。

第一類(lèi)是基于高斯消去法的直接求解器。直接計(jì)算方法可分為填充元優(yōu)化、符號(hào)分解、數(shù)值分解、消元回代4步。其中填充元優(yōu)化是對(duì)矩陣進(jìn)行重排序，在保證數(shù)值穩(wěn)定性的同時(shí)，使得后續(xù)矩陣分解時(shí)產(chǎn)生盡可能少的填充元，從而提高計(jì)算效率。常見(jiàn)的方法如最小度算法、重復(fù)剖分算法以及多重剖分算法。而符號(hào)分解則是現(xiàn)代快速直接解法特有的技術(shù)，通過(guò)使用消去樹(shù)這一數(shù)據(jù)結(jié)構(gòu)，確定分解因子矩陣的非零元位置，指導(dǎo)分解路徑，如基于一般消去樹(shù)的多波前方法和基于超節(jié)點(diǎn)消去樹(shù)的多元最小度方法等。直接解法的優(yōu)點(diǎn)是計(jì)算穩(wěn)定可靠，因此也成為大量通用有限元計(jì)算軟件中的默認(rèn)求解器；其缺點(diǎn)是計(jì)算過(guò)程中產(chǎn)生的分解因子、波前矩陣等中間變量規(guī)模不定，隨著問(wèn)題規(guī)模增大，超出內(nèi)存容量后超出部分只能存儲(chǔ)在硬盤(pán)中，硬盤(pán)的I/O操作成為計(jì)算瓶頸制約了直接解法的性能。

第二類(lèi)是迭代解法。預(yù)計(jì)于Gauss消去法的直接解法不同，迭代解法是在特定的解空間中不斷尋找更接近真解的元素，其中最常見(jiàn)的是基于krylov子空間的迭代方法：如先驗(yàn)共軛梯度迭代算法（preconditioned conjugate gradient，PCG）方法、雙共軛梯度（biconjugate gradient，BiCG）方法等。在迭代求解過(guò)程中需要不斷進(jìn)行稀疏矩陣向量乘法操作（sparse matrix-vector multiplication，SpMV）來(lái)獲得當(dāng)前解對(duì)應(yīng)的殘差，并基于此進(jìn)行修正。如何高效地執(zhí)行SpMV操作成為該類(lèi)解法的研究重點(diǎn)。研究人員通過(guò)調(diào)整稀疏矩陣的存儲(chǔ)格式來(lái)實(shí)現(xiàn)內(nèi)存的合并讀寫(xiě)，優(yōu)化SpMV操作的性能。常見(jiàn)的存儲(chǔ)格式可分為兩類(lèi)：行主元的COO格式、CSR格式、CTO格式以及列主元的ELL格式，Sliced ELL格式等［29］。針對(duì)硬件配置選擇特定的存儲(chǔ)格式，在存儲(chǔ)空間最小化同時(shí)通過(guò)緩存的合并讀寫(xiě)實(shí)現(xiàn)時(shí)空效率的最大化。

迭代求解器中另一個(gè)關(guān)鍵技術(shù)是預(yù)條件處理。隨著問(wèn)題規(guī)模增大，剛性矩陣的條件數(shù)越來(lái)越大，通過(guò)預(yù)條件處理能夠顯著提升迭代算法的收斂性能。在并行求解器中，剛性矩陣往往是分塊存儲(chǔ)，如何根據(jù)局部矩陣信息生成預(yù)條件矩陣是并行求解器研究的熱點(diǎn)，目前發(fā)展的并行預(yù)條件矩陣有：不完全LU分解預(yù)條件子，分塊Jacobi預(yù)條件子，最小二乘多項(xiàng)式預(yù)條件子以及多重網(wǎng)格預(yù)條件子等。相較于直接求解器，迭代求解器計(jì)算過(guò)程中對(duì)存儲(chǔ)的需求固定，借助區(qū)域分解方法容易實(shí)現(xiàn)稀疏矩陣-向量乘法操作的并行化，適合求解較大規(guī)模問(wèn)題；而迭代解法的缺點(diǎn)也非常明顯：無(wú)法保證收斂、生成性態(tài)優(yōu)良的預(yù)條件矩陣較為耗時(shí)。

縱觀有限元發(fā)展歷史，求解器經(jīng)歷了迭代方法（Gauss-Seidel方法）、直接方法（變帶寬方法）、PCG方法再到多波前等直接解法的螺旋式發(fā)展［30］。如今直接解法能夠有效求解百萬(wàn)自由度的有限元方程，而迭代解法多用于更大規(guī)模問(wèn)題［31］。無(wú)論是直接解法還是迭代解法，性能的提升都是通過(guò)發(fā)掘稀疏矩陣自身結(jié)構(gòu)并與計(jì)算機(jī)硬件相結(jié)合來(lái)實(shí)現(xiàn)，隨著計(jì)算機(jī)硬件發(fā)展，如何利用多核和眾核處理器實(shí)現(xiàn)高效求解成為新的問(wèn)題。

3 高性能有限元技術(shù)發(fā)展現(xiàn)狀與展望

3.1 從工業(yè)軟件角度重新認(rèn)識(shí)電磁仿真軟件

隨著物聯(lián)網(wǎng)（internet of things，IoT）技術(shù)的成熟，實(shí)時(shí)獲取設(shè)備運(yùn)行過(guò)程中的各種參數(shù)成為可能。為了有效利用這些參數(shù)對(duì)設(shè)備全生命周期進(jìn)行預(yù)測(cè)和控制，數(shù)字孿生（digital twin）這一概念被提出。由于實(shí)際設(shè)備的物理模型過(guò)于復(fù)雜，傳統(tǒng)數(shù)值仿真往往采用簡(jiǎn)化模型：如使用集總參數(shù)模型、將非線(xiàn)性時(shí)變的材料特性線(xiàn)性化等。仿真模型和實(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)行過(guò)程中出現(xiàn)的故障進(jìn)行檢測(cè)、定位以及預(yù)測(cè)。

數(shù)字孿生概念的出現(xiàn)極大地拓寬了電磁場(chǎng)數(shù)值仿真的應(yīng)用場(chǎng)景，同時(shí)也提出新的挑戰(zhàn)：精細(xì)的模型往往兼具時(shí)域、三維、非線(xiàn)性等特征，某些復(fù)雜的問(wèn)題還涉及熱場(chǎng)、流場(chǎng)、應(yīng)力場(chǎng)等多物理場(chǎng)耦合。如何利用邊緣計(jì)算、云計(jì)算等新興技術(shù)實(shí)現(xiàn)穩(wěn)定高效地求解是一個(gè)亟待解決的問(wèn)題。

3.2 借力國(guó)際開(kāi)放資源

在國(guó)際科學(xué)計(jì)算社區(qū)中不同研究團(tuán)隊(duì)分別提出并發(fā)布了用于偏微分方程求解的高性能數(shù)值計(jì)算庫(kù)：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專(zhuān)門(mén)用于求解偏微分方程；Trillinos使用面向?qū)ο蟮能浖蚣?，通過(guò)調(diào)用各個(gè)組件實(shí)現(xiàn)大規(guī)模復(fù)雜多物理場(chǎng)求解；另外還存在大量基于上述計(jì)算組件二次開(kāi)發(fā)的有限元計(jì)算軟件，如Deal.Ⅱ、FENICs、MFEM等。這些開(kāi)源的高性能計(jì)算軟件是當(dāng)前偏微分方程數(shù)值解研究的重要基礎(chǔ)，在計(jì)算流體力學(xué)、計(jì)算力學(xué)等領(lǐng)域都得到廣泛的應(yīng)用，但在計(jì)算電磁學(xué)領(lǐng)域未見(jiàn)報(bào)道。

3.3 結(jié)合人工智能技術(shù)的數(shù)值計(jì)算方法

近年來(lái)，大數(shù)據(jù)和人工智能技術(shù)突飛猛進(jìn)，深度神經(jīng)網(wǎng)絡(luò)技術(shù)給傳統(tǒng)科學(xué)計(jì)算帶來(lái)新的詮釋。神經(jīng)計(jì)算框架之下，科學(xué)計(jì)算問(wèn)題可分為3類(lèi)［32］：第1類(lèi)問(wèn)題是完全由數(shù)據(jù)驅(qū)動(dòng)的問(wèn)題，常見(jiàn)的如對(duì)風(fēng)、光等自然資源的實(shí)時(shí)預(yù)測(cè)。由于實(shí)際模型未知，往往直接簡(jiǎn)化為黑箱模型進(jìn)行處理，采用神經(jīng)網(wǎng)絡(luò)對(duì)數(shù)據(jù)進(jìn)行擬合、挖掘；第2類(lèi)問(wèn)題則是基于部分物理機(jī)理和少量測(cè)量數(shù)據(jù)的不確定性模型，常見(jiàn)的如逆散射問(wèn)題，即控制方程已知，根據(jù)測(cè)量數(shù)據(jù)反演出待求的模型參數(shù)。這一類(lèi)問(wèn)題往往可以寫(xiě)為PDE約束下的優(yōu)化問(wèn)題，借助神經(jīng)網(wǎng)絡(luò)這一優(yōu)化求解器進(jìn)行求解；第3類(lèi)問(wèn)題是基于已知模型的確定性問(wèn)題，經(jīng)典計(jì)算電磁學(xué)問(wèn)題往往屬于這一類(lèi)。有限元過(guò)程可以等價(jià)改寫(xiě)為神經(jīng)網(wǎng)絡(luò)的形式（FEA-net），即每個(gè)神經(jīng)元的輸入為單元幾何信息，輸出為單元?jiǎng)傂跃仃嚕?3］。應(yīng)用神經(jīng)網(wǎng)絡(luò)技術(shù)求解上述3類(lèi)問(wèn)題具有兩個(gè)顯著優(yōu)勢(shì)：一是借助各類(lèi)演化算法能夠處理復(fù)雜的非線(xiàn)性?xún)?yōu)化問(wèn)題；二是神經(jīng)網(wǎng)絡(luò)本身是對(duì)計(jì)算任務(wù)的細(xì)粒度分解，能夠?qū)崿F(xiàn)高性能計(jì)算。將經(jīng)典的計(jì)算電磁學(xué)問(wèn)題與如今蓬勃發(fā)展的神經(jīng)計(jì)算相結(jié)合，也是未來(lái)的發(fā)展方向之一。

3.4 發(fā)展異構(gòu)計(jì)算和跨平臺(tái)的高性能計(jì)算方法

高性能計(jì)算技術(shù)的發(fā)展與高性能計(jì)算機(jī)的發(fā)展是密切相關(guān)的，基于高性能計(jì)算技術(shù)的有限元方法在計(jì)算電磁學(xué)中的應(yīng)用經(jīng)歷了從最開(kāi)始只有單核CPU，發(fā)展到多核處理器，再到多核處理器和眾核處理器并存的異構(gòu)計(jì)算環(huán)境。隨著硬件平臺(tái)在不斷更新?lián)Q代，單位計(jì)算資源的成本越來(lái)越低：1997年，提供1TFlops計(jì)算能力的硬件需要花費(fèi)4 600萬(wàn)美元，而在2017年，使用相同的計(jì)算資源僅需要30美元［34］。相較于發(fā)展精巧的計(jì)算方法，高效利用“廉價(jià)”的硬件資源求解復(fù)雜的電磁學(xué)問(wèn)題更為重要，也更加有效。

在目前高性能有限元計(jì)算方案中，還存在以下問(wèn)題：其一，現(xiàn)有研究針對(duì)特定硬件的加速技術(shù)，這些技術(shù)往往受限于特定的問(wèn)題或特定的計(jì)算平臺(tái)，尚缺乏一款具有可移植性、可拓展性以及較高硬件利用率的計(jì)算程序；其二，除多核與眾核處理器外，張量處理器 TPU［35］、神經(jīng)計(jì)算單元 NPU［36］、量子處理器［37］等新興的計(jì)算單元也在不斷發(fā)展，有限元方法與這些硬件的結(jié)合也是一個(gè)亟待解決的問(wèn)題。

4 結(jié)語(yǔ)

有限元方法作為一種通用的數(shù)值求解技術(shù)，以其強(qiáng)大的計(jì)算能力，在電磁仿真領(lǐng)域獲得巨大成功。由于新應(yīng)用場(chǎng)景的出現(xiàn)，電力系統(tǒng)高壓設(shè)備電磁仿真對(duì)有限元方法提出了新的需求。從技術(shù)路徑發(fā)展上，高性能有限元技術(shù)是上述問(wèn)題的一個(gè)可能的解決方案。本文介紹了高性能有限元方法的主要技術(shù)要點(diǎn)，討論其相關(guān)的發(fā)展途徑，希望對(duì)高性能有限元方法在電力系統(tǒng)高壓設(shè)備仿真的應(yīng)用有所推動(dòng)。