肖調(diào)杰,周 峰,鄭翾宇,劉 劍,陳 琳,劉 杰,易明寬, 陳旭光,龔春葉,4,楊 博,甘新標(biāo),李勝國,左 克

(1.國防科技大學(xué)計(jì)算機(jī)學(xué)院,湖南 長沙 410073;2.國防科技大學(xué)高端裝備數(shù)字化軟件實(shí)驗(yàn)室,湖南 長沙 410073; 3.東華理工大學(xué)地球物理與測控技術(shù)學(xué)院,江西 南昌 330013;4.國家超級計(jì)算天津中心,天津 300457)

1 引言

地電磁學(xué)中的頻率域電磁法是一類重要的地球物理方法,其通過觀測由不同頻率天然場源或人工場源激發(fā)的電磁場在地下介質(zhì)中產(chǎn)生的感應(yīng)二次場,并分析其分布規(guī)律,進(jìn)而獲得地下電性分布特征及相關(guān)電性結(jié)構(gòu)參數(shù)[1-4]。頻率域電磁法種類繁多,包含大地電磁測深法、可控源音頻大地電磁法、海洋可控源電磁法等主動源方法以及大地電磁法、音頻大地電磁法、甚低頻法等被動源方法,具有探測深度大、分辨率高及抗干擾能力強(qiáng)等優(yōu)點(diǎn),在地球深部結(jié)構(gòu)探測、地?zé)峥碧?、礦產(chǎn)與油氣勘查及環(huán)境與工程勘探等諸多領(lǐng)域有著廣泛的應(yīng)用[5-10]。

數(shù)值方法和線性方程組求解是頻率域電磁法數(shù)值模擬的2個(gè)重要方面。一方面,有限單元法、有限差分法和積分方程法等是當(dāng)前頻率域電磁法中使用最為廣泛的3種數(shù)值模擬方法[11-17]。其中有限單元法和有限差分法屬于微分方程法,基于Maxwell方程的微分形式,求解時(shí)需要將全域作為計(jì)算區(qū)域,會造成較大的未知量和計(jì)算量,當(dāng)模型規(guī)模較大時(shí),其計(jì)算自由度可迅速增加至數(shù)千萬、數(shù)億量級甚至更大。積分方程法基于Maxwell方程的積分形式,是數(shù)值模擬算法中應(yīng)用得最早的一種方法。采用積分方程法計(jì)算電磁響應(yīng)時(shí),只需對目標(biāo)異常體進(jìn)行剖分,相比于有限單元法和有限差分法,積分方程所需求解的未知量大大減少。由于引入了格林函數(shù),積分方程法避免了微分方程法中場值傳遞造成的累積誤差,并且不需要施加吸收邊界條件或截?cái)噙吔鐥l件,從而避免了截?cái)噙吔缯`差,其數(shù)值結(jié)果具有半解析解的較高精度[18,19]。另一方面,相比于迭代法獲得近似解,直接解法可以獲得精確解。在三維數(shù)值模擬中,線性方程組求解占據(jù)了大部分計(jì)算量,求解速度直接關(guān)系到三維數(shù)值模擬的計(jì)算速度。目前方程組求解可以分為直接解法和迭代解法。迭代解法通過迭代過程中的不斷修正來獲得線性方程組的近似解。相比迭代解法,直接解法計(jì)算量非常大,且需要大量存儲空間,但是對于解存在的線性方程組,直接解法總能夠在有限的時(shí)間內(nèi)獲得其精確解,且具有非常高的穩(wěn)定性[20,21]。

當(dāng)前,頻率域電磁積分方程法數(shù)值模擬存在內(nèi)存受限和計(jì)算速度慢等問題,而大規(guī)模模型多頻點(diǎn)、多場點(diǎn)的計(jì)算已成為其實(shí)際應(yīng)用中的瓶頸。隨著Harrington[22]在1968年提出了矩量法,此后積分方程法在計(jì)算電磁領(lǐng)域得到了快速發(fā)展。而在三維地電領(lǐng)域,Hohmann[23]在1975年基于六面體網(wǎng)格闡述了奇異核積分技術(shù),開啟了將積分方程法用于求解三維地電問題的大門,隨后Weidelt[24]推導(dǎo)得到了層狀介質(zhì)張量格林函數(shù)的積分表達(dá)式,奠定了積分方程法在地電問題中的應(yīng)用基礎(chǔ)。之后,諸多研究人員對積分方程法在頻率域電磁法中的應(yīng)用開展了大量的研究[10,12]。積分方程法在求解過程中需要求解稠密線性方程組,若求解區(qū)域的未知量為W,則其存儲需量為O(W2),直接求解的計(jì)算復(fù)雜度為O(W3),隨著計(jì)算規(guī)模的增大,面臨著內(nèi)存受限和計(jì)算速度慢等問題。在優(yōu)化存儲方面,Ting等[25]提出利用格林函數(shù)的對稱性來減少系數(shù)矩陣的存儲,但是隨著模型網(wǎng)格數(shù)量的進(jìn)一步增大,依然面臨著內(nèi)存受限的問題。在提高計(jì)算速度方面,Wannamaker等[26]基于大網(wǎng)格單元的格林系數(shù)插值出小網(wǎng)格單元的格林系數(shù)來加快獲得格林系數(shù)矩陣;Aksun[27]、袁建生等[28]采用復(fù)鏡像法加快格林函數(shù)的計(jì)算。為加快方程組的求解速度,Tripp等[29]把系數(shù)矩陣拆分成了若干個(gè)小矩陣,Millard等[30]提出了穩(wěn)定雙共軛梯度的快速傅里葉變換法。此外,Zhadnov等[31,32]、Abubakar等[33]、Ueda等[34]和Kruglyakov等[35]深入研究了Born近似、擬線性近似等方法,以降低求解精度為代價(jià)提高了計(jì)算速度。而隨著計(jì)算規(guī)模、場點(diǎn)和頻點(diǎn)的迅速增加,所需內(nèi)存和計(jì)算量將急劇遞增。上述研究在一定程度上緩解了內(nèi)存受限和計(jì)算速度慢的問題,但未能從根本上給出合理的解決方案。

因此,本文面向頻率域電磁場數(shù)值模擬,采用積分方程法、分布式存儲和直接解法等技術(shù),設(shè)計(jì)了一種多層級多粒度混合并行策略,解決了內(nèi)存受限和計(jì)算速度慢等問題,實(shí)現(xiàn)了頻點(diǎn)間并行、阻抗矩陣并行填充、方程組并行求解的快速、高精度、高可擴(kuò)展性的頻率域電磁三維積分方程數(shù)值模擬方法。

2 物理模型

圖1所示為頻率域電磁法三維地電模型示意圖。在地表有一接地導(dǎo)線作為電性發(fā)射源(Transmitter),源的長度為L,發(fā)射電流為I。積分方程法求解范圍是電磁參數(shù)區(qū)別于背景介質(zhì)的異常區(qū)域,圖1中所示背景區(qū)域是電導(dǎo)率為σb的均勻半空間,其中存在電導(dǎo)率為σ的異常區(qū)域。給定發(fā)射源、介質(zhì)電性參數(shù)及空間幾何等信息,可通過數(shù)值模擬計(jì)算得到整個(gè)空間的電磁場分布。

Figure 1 Diagram of geo-electromagnetic model圖1 三維地電模型示意圖

3 積分方程法

研究任何形式的電磁場都必須從Maxwell方程組出發(fā),電磁場通過傅里葉變換,可將時(shí)間域的電磁場轉(zhuǎn)換為一系列諧變場的組合,取時(shí)諧因子為e-iwt,在準(zhǔn)靜態(tài)條件下,Maxwell方程組變?yōu)槭?1)～式(4)[5]:

×E=iωμΗ

(1)

×H=(σ-iωε)E+Js

(2)

·εE=ρ

(3)

·μH=0

(4)

其中,E是電場強(qiáng)度矢量,H是磁場強(qiáng)度矢量,ρ是自由電荷密度,ε是介質(zhì)的介電常數(shù),μ是介質(zhì)磁導(dǎo)率,σ是介質(zhì)的電導(dǎo)率,Js是激發(fā)源,ω=2πf是角頻率,f是頻率,i是虛數(shù)單位。

對式(1)求旋度后,結(jié)合式(2),可得到電場雙旋度方程,如式(5)所示:

××E-iωμ(σ-iωε)E=iωμJs

(5)

根據(jù)疊加原理,將總電場分解為入射電場Einc和散射電場Esca,有:

E=Einc+Esca

(6)

采用電張量格林函數(shù),并結(jié)合式(5),可以得到積分公式,如式(7)所示:

(7)

其中,Ωs是場源積分區(qū)域,Ω是異常體積分區(qū)域,G是并矢Green函數(shù)。當(dāng)背景模型為均勻半空間時(shí),G是均勻半空間并矢Green函數(shù);當(dāng)背景模型為層狀介質(zhì)時(shí),G是層狀介質(zhì)并矢函數(shù)[25,26]。

在進(jìn)行數(shù)值模擬之前,需要對式(7)進(jìn)行剖分離散。本文采用四面體單元,將異常體離散成N個(gè)四面體子單元。采用適用于介質(zhì)的體積分方程法,假設(shè)每個(gè)單元內(nèi)部的電導(dǎo)率和電場保持不變,取為單元中心ri(i=1,2,…,M)處的值(如圖2所示,其中Ex,Ey和Ez分別為x,y和z方向的電場強(qiáng)度),可得式(8):

E(rm)=Einc(rm)+

(8)

其中,Ei為單元中心處的電場。

Figure 2 Diagram of tetrahedral element圖2 四面體單元示意圖

將式(8)寫成矩陣形式,可形成關(guān)于異常體區(qū)域3N個(gè)未知電場的線性方程組,如式(9)～式(12)所示:

A·E=Einc

(9)

(10)

(11)

(12)

求解式(8),可得到異常體區(qū)域的電場值。而最終某一點(diǎn)的總場由該點(diǎn)的入射場和散射場疊加而成,再通過式(8)即可獲取空間中任一點(diǎn)處的電場值。

4 并行策略

4.1 多級并行架構(gòu)

在計(jì)算之前,采用國產(chǎn)網(wǎng)格剖分軟件YH-Grid或開源軟件Tetgen對模型進(jìn)行網(wǎng)格離散,并準(zhǔn)備好相應(yīng)的輸入文件。電磁法的頻點(diǎn)間計(jì)算無數(shù)據(jù)依賴關(guān)系,具有天然的可并行性。對單頻點(diǎn)的串行算法進(jìn)行計(jì)算時(shí)間占比分析,發(fā)現(xiàn)其阻抗矩陣填充、方程組求解和觀測點(diǎn)場值計(jì)算占總計(jì)算時(shí)間90%以上。圖3為本文設(shè)計(jì)的并行算法流程示意圖。首先,采取進(jìn)程分組的方式將不同頻點(diǎn)劃分到不同的進(jìn)程組。然后,各進(jìn)程組的主進(jìn)程讀取輸入文件并廣播給組內(nèi)所有進(jìn)程。接著,并行填充阻抗矩陣A,采用分布式存儲。最后加入邊界條件,對方程組進(jìn)行并行求解。求解方程組之后,可進(jìn)一步對觀測點(diǎn)的場值進(jìn)行并行計(jì)算。最后進(jìn)行數(shù)據(jù)收集和后處理。

Figure 3 Flow diagram of multi-level parallel algorithm圖3 多級并行算法流程示意圖

4.2 負(fù)載均衡

為保持良好的負(fù)載均衡與可擴(kuò)展性,采取分布式存儲,采用廣泛應(yīng)用在ScaLAPACK、SLATE、MAGMA等稠密線性代數(shù)庫中的二維Block- cyclic數(shù)據(jù)映射方法。假設(shè)一共有N個(gè)四面體單元,則阻抗矩陣A的大小為3N×3N。將矩陣元素以塊的形式不重復(fù)地循環(huán)分配/存儲在P×Q個(gè)不同進(jìn)程之中,其中P和Q分別為二維進(jìn)程網(wǎng)格中垂直方向和水平方向的進(jìn)程個(gè)數(shù)。二維進(jìn)程網(wǎng)格如圖4所示,在邏輯上將9N2個(gè)矩陣元素分成一個(gè)個(gè)大小為mb×nb的分塊小矩陣,其中每個(gè)格子代表一個(gè)分塊小矩陣,相同顏色的分塊小矩陣存儲在同一個(gè)進(jìn)程中,每個(gè)進(jìn)程保存的矩陣元素個(gè)數(shù)為9N2/(P×Q),一般要求緩存M≤9N2/(P×Q)。

Figure 4 Block-cyclic mapping圖4 Block-cyclic映射

設(shè)當(dāng)前進(jìn)程在進(jìn)程網(wǎng)格中的索引為(my_row,my_col),圖4a中矩陣元素的全局索引為(ig,jg),圖4b中當(dāng)前進(jìn)程中矩陣元素的局部索引為(il,jl),索引從0開始,則有:

my_row=(ig/mb)%P

(13)

my_col=(jg/nb)%Q

(14)

當(dāng)前進(jìn)程局部索引和全局索引之間的關(guān)系滿足式(15)～式(18):

il=(ig/mb)/P×mb+ig%mb

(15)

jl=(jg/nb)/Q×nb+jg/nb

(16)

ig=P×mb×(il/mb)+il%mb

(17)

jg=Q×nb×(jl/nb)+jl%nb

(18)

4.3 網(wǎng)格分發(fā)

根據(jù)式(10)可知,阻抗矩陣A由四面體兩兩之間的相互關(guān)系系數(shù)矩陣組合而成,如圖5所示,每個(gè)四面體單元都對應(yīng)3行和3列。

Figure 5 Impedance matrix A圖5 阻抗矩陣A

以總進(jìn)程為6、P×Q=2×3、阻抗矩陣A分成9×9個(gè)分塊小矩陣Tij(i=0,…,8;j=0,…,8)為例,每個(gè)分塊小矩陣大小為NB×NB,Tij如式(19)所示:

(19)

如圖6所示,按照塊循環(huán)(block-cyclic)方式將阻抗矩陣A映射到2×3的進(jìn)程網(wǎng)格中。

Figure 6 Impedance matrix圖6 阻抗矩陣分布示意圖

如圖7所示,將每行或每列塊(block)對應(yīng)的單元網(wǎng)格信息存放在一個(gè)包(bucket)中,矩陣塊是方陣,第i行block和第i列block涉及到的網(wǎng)格單元相同,對應(yīng)同一個(gè)bucket,在圖中不做展示。其中每個(gè)四面體需要保存其全局編號(index)、屬性編號(id)、4個(gè)節(jié)點(diǎn)位置信息(P0(x0,y0,z0),P1(x1,y1,z1),P2(x2,y2,z2),P3(x3,y3,z3)),即2個(gè)整數(shù)和12個(gè)浮點(diǎn)數(shù),分別為(index,id,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3)。

Figure 7 Grid distribution圖7 網(wǎng)格分發(fā)

根據(jù)圖6b,進(jìn)程(0,0)存儲的矩陣塊如圖8所示,進(jìn)程(0,0)所存儲的矩陣塊為T0,0,T0,3,T0,6,T2,0,T2,3,T2,6,T4,0,T4,3,T4,6,T6,0,T6,3,T6,6,T8,0,T8,3,T8,6,填充Tij需要第i和第j個(gè)bucket,其中行涉及到的bucket為第0,2,4,6和8個(gè)bucket,列涉及到的bucket為第0,3和6個(gè)bucket。因此,一起需要第0,2,3,4,6和8個(gè)bucket。將行和列涉及到的bucket指針分別保存在rowGrid和colGrid2個(gè)不同的數(shù)組中。一般地,第i個(gè)bucket由進(jìn)程(i%P,i%Q)進(jìn)行廣播。首先,負(fù)責(zé)讀入網(wǎng)格和切分網(wǎng)格到bucket的主進(jìn)程將第i個(gè)bucket發(fā)送給進(jìn)程(i%P,i%Q),再由該進(jìn)程在row_comm和col_comm分別廣播,行進(jìn)程號等于i%P的進(jìn)程將收到的bucket指針保存在rowGrid中。列進(jìn)程號等于i%Q的進(jìn)程將收到的bucket指針保存在colGrid中,即完成網(wǎng)格分發(fā)。

其中,負(fù)責(zé)廣播的進(jìn)程保存的bucket指針同時(shí)保存在rowGrid和colGrid中。本文引入isBoth的哈希表進(jìn)行判斷,以避免內(nèi)存重復(fù)釋放。

Figure 8 Matrix block of process (0,0)圖8 進(jìn)程(0,0)存儲的矩陣塊

4.4 阻抗矩陣填充

網(wǎng)格分發(fā)完成之后,各進(jìn)程已得到所需網(wǎng)格信息。接下來,只需要在rowGrid和colGrid中對每個(gè)四面體分別進(jìn)行計(jì)算,然后保存到本地阻抗矩陣即可。假設(shè)rowGrid有m個(gè)bucket,colGrid有n個(gè)bucket,阻抗矩陣填充如算法1所示[36]。

算法1阻抗矩陣填充

輸入：m個(gè)bucket組成的網(wǎng)格塊rowGrid和n個(gè)bucket組成的網(wǎng)格塊colGrid。

輸出：本地阻抗矩陣A。

forj=0:ndo

fori=0:mdo

Aptr=A+i*NB+j*NB*ld;/*跳轉(zhuǎn)到對應(yīng)NB塊的首地址*/

nElement←colGrid的四面體個(gè)數(shù);

mElement←rowGrid的四面體個(gè)數(shù);

fork=0:nElementdo

計(jì)算本地填充的列方向開始位置n0和結(jié)束位置n1;

forp=0:mElementdo

計(jì)算本地填充的行方向開始位置m0和結(jié)束位置m1;

計(jì)算四面體之間產(chǎn)生的矩陣T;

Aptr(m0:m1,n0:n1)←

T(m0:m1,n0:n1);

endfor

endfor

endfor

endfor

4.5 復(fù)數(shù)稠密線性系統(tǒng)LU分解

積分方程法最終形成一個(gè)線性稠密的復(fù)數(shù)線性方程組,根據(jù)式(10)可知,采用部分選主元的LU分解進(jìn)行求解[37],如式(20)所示:

(20)

其中,A11∈CNB×NB,A12∈CNB×(N-NB),A21∈C(N-NB)×NB,A22∈C(N-NB)×(N-NB)。

4.6 觀測場值并行計(jì)算

已知異常體各離散單元的場值,便可根據(jù)式(8)求得任意觀測點(diǎn)處的場值。各觀測點(diǎn)的場值為背景場與二次場之和,其中二次場為地下異常體在該處產(chǎn)生的感應(yīng)場值,由異常體中各四面體在該處產(chǎn)生的場值累加組成。當(dāng)模型規(guī)模較大、觀測點(diǎn)數(shù)較多時(shí),觀測點(diǎn)處的場值計(jì)算耗時(shí)巨大?？紤]到實(shí)際應(yīng)用場景,進(jìn)程數(shù)通常少于觀測點(diǎn)數(shù),因此在各進(jìn)程組/單個(gè)頻點(diǎn)內(nèi)對觀測點(diǎn)處的場值進(jìn)行并行計(jì)算。采用主從模式的并行加速策略,各進(jìn)程組內(nèi)主進(jìn)程負(fù)責(zé)任務(wù)調(diào)度和結(jié)果收集,從進(jìn)程負(fù)責(zé)計(jì)算并向主進(jìn)程返回計(jì)算結(jié)果。計(jì)算流程如圖9所示。

Figure 9 Flow chart of parallel calculation of field at observation points圖9 觀測點(diǎn)場值并行計(jì)算流程圖

5 驗(yàn)證與測試

5.1 正確性和準(zhǔn)確度驗(yàn)證

如圖10所示,在電導(dǎo)率為0.02 S/m的均勻半空間中有一個(gè)電導(dǎo)率為0.2 S/m的低阻體。低阻體的埋深為100 m,尺寸為120 m × 200 m × 400 m,中心坐標(biāo)為(0 m,1 000 m,300 m)。偶極子源沿y方向,長度為1 m,發(fā)射電流為100 A,頻率為3 Hz,中心坐標(biāo)為(0 m,0 m,0 m)。沿y軸方向在異常體上方布設(shè)間距為20 m的41個(gè)觀測點(diǎn),起點(diǎn)位置為(0 m,600 m,0 m),終點(diǎn)位置為(0 m,1 400 m,0 m)。

Figure 10 Testing model:There is a conductive body embedded in a half-space圖10 測試模型:均勻半空間中存在一低阻體

選取Ey的實(shí)部、虛部與前人結(jié)果[38]進(jìn)行對比,如圖11所示,可以看到,實(shí)部和虛部的相對誤差在6%以內(nèi)。

Figure 11 Comparison of electric field amplitudes圖11 電場振幅對比

Figure 13 Distribution of real component of Ey圖13 Ey實(shí)部分布

5.2 性能測試

以可控源為例,如圖12所示。按對數(shù)在10-1～102Hz均勻取16個(gè)頻點(diǎn),共16 × 12495個(gè)未知量。若采用串行算法中的全存儲策略而不是本文中的分布式存儲,阻抗矩陣所需存儲空間大小為16×(12495×2)2×8/10243=74 GB,而隨著問題規(guī)模的增加,所需存儲空間會進(jìn)一步極速劇增。因此,采用本文算法,在天河高性能計(jì)算系統(tǒng)上進(jìn)行測試,每個(gè)節(jié)點(diǎn)分配32個(gè)進(jìn)程,共設(shè)置861個(gè)觀測點(diǎn)。測試了1個(gè)節(jié)點(diǎn)32個(gè)進(jìn)程到256個(gè)節(jié)點(diǎn)8 192個(gè)進(jìn)程的不同規(guī)模。

Figure 12 Large-scale model for testing圖12 大規(guī)模測試模型

參考Dublin Test Model,設(shè)計(jì)組合體模型如圖12所示。在圖12中,電導(dǎo)率σb為0.001 25 S/m的均勻半空間存在一個(gè)埋深500 m的低阻異常組合體,組合體的電導(dǎo)率σ1、σ2及σ3分別為0.02 S/m,0.01 S/m和0.005 S/m,對應(yīng)塊體的尺寸分別為300 m×400 m×600 m,300 m×400 m×600 m,2000 m×400 m×400 m,發(fā)射場源為沿y方向布設(shè),長度100 m,電流大小為2 A,觀測區(qū)域位于異常體的正上方,x方向:-2 000 m～2 000 m,y方向:-1 000 m～1 000 m,x和y方向間距100 m,共布設(shè)861個(gè)觀測點(diǎn)。坐標(biāo)系的源點(diǎn)O設(shè)定在異常組合體在水平面投影的中心點(diǎn)處。

Figure 14 Distribution of imaginary component of Ey圖14 Ey虛部分布

計(jì)算結(jié)果如圖13和圖14所示,分別展示了Ey的實(shí)部和虛部,可以看到異常體上方的電場有明顯的變化,不同頻點(diǎn)對異常體均有所反映。

可擴(kuò)展性測試結(jié)果如表1和圖15所示?？梢钥吹?本文開發(fā)的并行程序具有較好的可擴(kuò)展性,相比1個(gè)節(jié)點(diǎn)32進(jìn)程,當(dāng)測試規(guī)模達(dá)到256節(jié)點(diǎn)8 192進(jìn)程時(shí),加速比為69.69,并行效率為27.22%。

Table 1 Speedup and parallel efficiency of large-scale model表1 大規(guī)模模型的加速比及并行效率

Figure 15 Speedup and parallel efficiency of large-scale model圖15 大規(guī)模模型的加速比和并行效率

6 結(jié)束語

本文詳細(xì)描述了一種頻率域電磁積分方程法大規(guī)模數(shù)值模擬程序的設(shè)計(jì)、實(shí)現(xiàn)方案和測試結(jié)果。通過對程序的驗(yàn)證和測試可以看出,該程序具有較高的精度和較好的可擴(kuò)展性,具備在大型超算平臺上大規(guī)模運(yùn)行的能力,可以滿足地電磁學(xué)中大地電磁、可控源音頻大地電磁等精細(xì)數(shù)值模擬的需求。