馬　駒， 肖調(diào)杰*， 王　赟

(1. 中國(guó)科學(xué)院　地球化學(xué)研究所，貴陽(yáng)　550002；2.中國(guó)科學(xué)院大學(xué)，北京　100049；3. 中國(guó)科學(xué)院　地質(zhì)與地球物理研究所，北京　100029)

?

MPI并行的節(jié)點(diǎn)大地電磁三維有限元正演

馬駒1,2， 肖調(diào)杰1,2*， 王赟1,3

(1. 中國(guó)科學(xué)院地球化學(xué)研究所，貴陽(yáng)550002；2.中國(guó)科學(xué)院大學(xué)，北京100049；3. 中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京100029)

摘要：用fortran語(yǔ)言編程實(shí)現(xiàn)了有限元三維大地電磁正演，通過(guò)層狀介質(zhì)模型、二維棱柱體模型及三維低阻體模型結(jié)果的對(duì)比，驗(yàn)證了所編寫(xiě)程序的正確性。首先通過(guò)加入第一類邊界條件，減少了最終求解方程組的維數(shù)，同時(shí)對(duì)系數(shù)矩陣的存儲(chǔ)采用非零存儲(chǔ)技術(shù)，大大降低了對(duì)計(jì)算機(jī)內(nèi)存的需求；最后在串行程序的基礎(chǔ)上，基于MPI實(shí)現(xiàn)了頻點(diǎn)間的并行，并對(duì)一個(gè)三維模型進(jìn)行計(jì)算，并行后開(kāi)啟4進(jìn)程時(shí)加速比達(dá)到了2.328，有效地減少了所需時(shí)間。

關(guān)鍵詞：大地電磁； 三維； 正演； 并行； MPI

0引言

自上世紀(jì)70年代以來(lái)很多學(xué)者都對(duì)電磁三維正演進(jìn)行了研究，主要有積分方程法[1-3]、有限差分法[4-6]及有限單元法[7]。Wannamaker[2]于對(duì)積分方程法的應(yīng)用于三維大地電磁正演的進(jìn)展進(jìn)行了論述；Mackie等[4]討論了有限差分法在三維大地電磁模擬中的應(yīng)用，并與前人的積分方程法結(jié)果進(jìn)行了對(duì)比。相對(duì)而言，有限元法更適用于地形復(fù)雜情況的三維電磁模擬， Thohru Mogi[7]利用二次場(chǎng)實(shí)現(xiàn)了三維正演，避免了電場(chǎng)在不連續(xù)界面不連續(xù)的問(wèn)題，并通過(guò)模型驗(yàn)證了其正確性； Yuji Mitsuhata等[8]實(shí)現(xiàn)了電矢量位磁標(biāo)量位三維有限元正演，效果不錯(cuò)； Puzyrev等[9]給出了基于磁矢量勢(shì)和電標(biāo)量勢(shì)的并行頻率域節(jié)點(diǎn)有限元法；Silva等[10]給出了基于多波前法的三維矢量有限元數(shù)值模擬結(jié)果。

目前，在地球物理學(xué)中，并行計(jì)算已經(jīng)得到越來(lái)越多的關(guān)注和應(yīng)用，當(dāng)前主要有兩種并行編程模式[11]，①基于消息傳遞型的進(jìn)程級(jí)并行編程模式MPI(Message Passing Interface)；②共享內(nèi)存型的線程級(jí)并行編程模式OpenMP(Open Multi-Processing)。相對(duì)于OpenMP ，這里選取移植性更強(qiáng)的MPI并行模式。作者基于有限單元法用fortran語(yǔ)言編程實(shí)現(xiàn)了大地電磁三維正演，采用六面體剖分，單元內(nèi)進(jìn)行線性插值，最后形成的總體剛度矩陣是一大型稀疏、帶狀、對(duì)稱但不正定的復(fù)系數(shù)矩陣[12]，最后形成的系數(shù)矩陣采用CSR(Compressed Sparse Row format)非零存儲(chǔ)以節(jié)省內(nèi)存空間，解方程采用高效、收斂速度快的雙共軛梯度法，同時(shí)采用直接法加入第一類邊界條件，降低了系數(shù)矩陣的維數(shù)[13]；最后在串行程序基礎(chǔ)上，基于MPI實(shí)現(xiàn)了頻點(diǎn)間的并行計(jì)算，有效地節(jié)省了時(shí)間。

1三維正演基本理論

1.1邊值問(wèn)題

取時(shí)間因子為e-iωt，則由Maxwell方程組可得式(1)。

(1)

其中：E是電場(chǎng)；k2=iωμσ+ω2εμ；ω為角頻率；μ是介質(zhì)的導(dǎo)磁系數(shù)；σ是介質(zhì)的電導(dǎo)率；ε是介質(zhì)的介電常數(shù)。

圖1　三維大地電磁正演模型示意圖Fig．1　Model of 3D magnetotelluric sounding

1.2六面體剖分

母單元是一個(gè)邊長(zhǎng)為2的正方體，取8個(gè)角點(diǎn)為節(jié)點(diǎn)，編號(hào)及坐標(biāo)如圖2(a)所示，圖2(b)為子單元。

圖2　六面體剖分單元Fig．2　Division of hexahedron(a)母單元；(b)子單元

構(gòu)造形函數(shù)為式(2)。

(2)

其中：ξi、ηi、γi是節(jié)點(diǎn)i在母單元中的坐標(biāo)。

1.3總體系數(shù)矩陣

由廣義變分原理建立泛函，如式(3)所示。

k2E·E]dV

(3)

加入散度條件：

1)空氣中散度條件為▽·E=0。

2)介質(zhì)中散度條件為▽·σE=0，得到式(4)。

在式(4)中加入邊界條件，進(jìn)行離散，然后在單元內(nèi)插值，再將積分轉(zhuǎn)換為線性函數(shù)，最后得到總體剛度矩陣，進(jìn)而得到線性方程組：

KE=0

(5)

解方程即可得各電場(chǎng)分量，進(jìn)而可根據(jù)電場(chǎng)分量計(jì)算得到各磁場(chǎng)分量。

1.4張量阻抗計(jì)算

(6)

(7)

(8)

(9)

式(7)定義的響應(yīng)稱為XY模式響應(yīng)，式(8)定義的響應(yīng)稱為YX模式響應(yīng)。進(jìn)一步便可計(jì)算得到相應(yīng)的視電阻率和相位如式(10)所示。

(10)

(i=x,y;j=x,y)

2模型計(jì)算

分別對(duì)層狀介質(zhì)模型、二維模型及三維模型進(jìn)行了計(jì)算并對(duì)比。

2.1三層層狀狀介質(zhì)模型

三層層狀介質(zhì)模型如圖3所示。我們用解析解與三維數(shù)值解進(jìn)行了對(duì)比，計(jì)算頻點(diǎn)為1 000 Hz、100 Hz、10 Hz、1 Hz、0.1 Hz、0.01 Hz與0.001 Hz共7個(gè)頻點(diǎn)。視電阻率曲線見(jiàn)圖4，相位曲線見(jiàn)圖5。其中視電阻率最大相對(duì)誤差為0.282 8%，相位最大相對(duì)誤差為0.048 23%，表明三維數(shù)值模擬結(jié)果與解析解擬合得非常好。

圖3　三層層狀介質(zhì)模型Fig．3　The model of layered media

圖4　三層層狀介質(zhì)三維大地電磁正演視電阻率值與解析解的對(duì)比Fig．4　The comparison of apparent resistivity between analytical results and numerical results

圖5　三層層狀介質(zhì)三維大地電磁正演相位與解析解的對(duì)比Fig．5　The comparison of phase between analytical results and numerical results

2.2二維棱柱體模型

二維棱柱體模型如圖6所示，棱柱體走向?yàn)閅方向，電阻率為10 Ω·m，幾何尺寸為3 000 m×6 000 m ，頂面埋深為3 000 m ，圍巖電阻率為100 Ω·m。

我們分別用二維有限單元法和三維有限單元法對(duì)該模型進(jìn)行了計(jì)算，計(jì)算頻點(diǎn)為3 Hz、1 Hz及0.5 Hz三個(gè)頻點(diǎn)。將三維YX模式、XY模式分別與二維TE模式、TM模式進(jìn)行了對(duì)比，圖7為三維YX模式與二維TE模式對(duì)比，二者結(jié)果擬合得非常好，最大相對(duì)誤差為0.378 7%，圖8為三維XY模式與二維TM模式對(duì)比，二者結(jié)果擬合得很好，最大相對(duì)誤差為1.288 8%。

圖6　低阻體模型Fig．6　Geometry of 2D prism model

圖7　YX模式與TE模式對(duì)比Fig．7　The comparison of YX mode and TE mode

圖8　XY模式與TM模式對(duì)比Fig．8　The comparison of XY mode and TM mode

2.3三維棱柱體

三維模型采用文獻(xiàn)[6]的三維棱柱體模型，如圖9所示，三維棱柱體尺寸為6 000 m×6 000 m×3 000 m，電阻率為10 Ω·m，頂面埋深為3 000 m，圍巖電阻率為100 Ω·m。計(jì)算頻點(diǎn)為3 Hz、1 Hz及0.1 Hz共3個(gè)頻點(diǎn)，將YX模式、XY模式分別與參考文獻(xiàn)[6]中的數(shù)值模擬結(jié)果進(jìn)行了對(duì)比，YX模式對(duì)比如圖10所示，可以看到不管是XY模式還是YX模式，3 Hz及1 Hz時(shí)二者擬合很好， 0.1 Hz時(shí)曲線形態(tài)基本一致，值有所差別，但最大相對(duì)誤差都在5%以內(nèi)。

圖9　三維低阻體模型Fig．9　Geometry of 3D prism model

圖10　YX模式視電阻率對(duì)比Fig．10　The apparent resistivity of YX mode

圖11　XY模式視電阻率對(duì)比Fig．11　The apparent resistivity of XY mode

通過(guò)以上三個(gè)模型的計(jì)算與對(duì)比，表明三維有限元數(shù)值模擬是準(zhǔn)確可靠的。

3MPI并行計(jì)算

3.1MPI并行實(shí)現(xiàn)

MPI(Massage Passing Interface)是目前最為通用的并行編程方式之一，也是分布式并行系統(tǒng)的主要編程環(huán)境，作者在實(shí)現(xiàn)了大地電磁三維正演基礎(chǔ)上，基于MPI實(shí)現(xiàn)了頻點(diǎn)間的并行，MPI并行算法偽代碼如下：

PROGRAM MT3DForwardFEM

USE MPI !調(diào)用MPI模塊

……

!初始化MPI系統(tǒng)；得到所有參加運(yùn)算的進(jìn)程的個(gè)數(shù)，放NUMPROCS中；得到當(dāng)前正在運(yùn)行的進(jìn)程的標(biāo)識(shí)號(hào)，放在MY_RANK中

CALL MPI_INIT(IERR)

CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NUMPROCS,IERR)

CALL MPI_COMM_RANK(MPI_COMM_WORLD,MY_RANK,IERR)

……

!為每個(gè)進(jìn)程都分配長(zhǎng)度一樣的一維數(shù)組用來(lái)存儲(chǔ)每個(gè)進(jìn)程得到的結(jié)果，F(xiàn)N是頻點(diǎn)數(shù)，!UNX是測(cè)點(diǎn)數(shù)，在這里統(tǒng)一分配(FN/NUMPROCS+1)*UNX大小，后面對(duì)結(jié)果再進(jìn)行一下處理，收集結(jié)果時(shí)就可以用MPI_GATHER()函數(shù)而不必使用麻煩的MPI_GATHERV()函數(shù)

ALLOCATE(RESXY((FN/NUMPROCS+1)*UNX))

ALLOCATE(APHXY((FN/NUMPROCS+1)*UNX))

ALLOCATE(RESYX((FN/NUMPROCS+1)*UNX))

ALLOCATE(APHYX((FN/NUMPROCS+1)*UNX))

……

IF(MY_RANK==0)THEN

讀取電阻率、坐標(biāo)等信息

END IF

CALL MPI_BCAST()!廣播主進(jìn)程讀取的輸入信息

……

DO

FI=MY_RANK+1,FN,NUMPROCS

… …

調(diào)用子程序計(jì)算形成總體系數(shù)矩陣、加入邊界條件值、解方程、計(jì)算視電阻率及相位

… …

END DO

……

!進(jìn)程同步后收集各進(jìn)程上的結(jié)果

CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)

!在主進(jìn)程上收集XY模式及YX模式的視電阻率及相位

CALL MPI_GATHER()

CALL MPI_GATHER()

CALL MPI_GATHER()

CALL MPI_GATHER()

……

寫(xiě)入結(jié)果

……

CALL MPI_FINALIZE(IERR) !退出MPI系統(tǒng)

END PROGRAM MT3DFowardFEM

程序在采用MPI_GATHER()函數(shù)的情況下實(shí)現(xiàn)了對(duì)任意頻點(diǎn)數(shù)的并行化。

3.2并行效率提升評(píng)價(jià)

對(duì)一個(gè)低阻體模型進(jìn)行了并行計(jì)算，模型如圖12所示，低阻體尺寸為800 m×800 m×1 050 m，頂面埋深為350 m，電阻率為40 Ω·m，圍巖電阻率為100 Ω·m。并行程序可以不加修改的在其他并行平臺(tái)上運(yùn)算，我們?cè)谒暮藱C(jī)上進(jìn)行運(yùn)算，配置為Intel(R) CoreTMi7-4790 CPU，主頻為3.6 GHz，內(nèi)存16 G，64位操作系統(tǒng)。從0.01 Hz～1 000 Hz范圍內(nèi)按對(duì)數(shù)均勻取26個(gè)頻點(diǎn)進(jìn)行并行計(jì)算，網(wǎng)格大小為42×42×23。并行計(jì)算時(shí)間見(jiàn)表1，圖13與圖14分別反映了開(kāi)啟不同進(jìn)程時(shí)加速比及效率的變化，其中，并行加速比=串行時(shí)間/并行計(jì)算時(shí)間；并行效率=并行加速比/并行進(jìn)程個(gè)數(shù)。

對(duì)圖13、圖14及表1進(jìn)行分析：當(dāng)開(kāi)啟進(jìn)程數(shù)增加時(shí)，并行效率降低；隨著進(jìn)程數(shù)量的增加，并行加速比繼續(xù)增大，但并行效率反而下降，這是因?yàn)殡S著進(jìn)程數(shù)的增加，用于進(jìn)程間通信的時(shí)間所占比例逐漸增大；1個(gè)進(jìn)程的并行計(jì)算時(shí)間比串行時(shí)間長(zhǎng)，這是因?yàn)椴捎貌⑿行枰~外的開(kāi)銷，總的看來(lái)，并行效果還是很明顯的，當(dāng)開(kāi)啟4個(gè)進(jìn)程時(shí)達(dá)到了2.3的加速比。

圖12　三維低阻體模型Fig．12　Geometry of 3D model

圖13　開(kāi)啟進(jìn)程數(shù)量不同時(shí)的加速比Fig．13　The speed-up with different 　numbers of threads

圖14　開(kāi)啟進(jìn)程數(shù)量不同時(shí)的效率Fig．14　The efficiency with different 　number of threads

3.3MPI并行計(jì)算結(jié)果

由于并行計(jì)算并沒(méi)有改變算法，因此并行結(jié)果和串行結(jié)果一樣。取X=-1 400 m、-800 m、-400 m、-200 m、0 m、200 m、400 m、800 m、1 400 m 測(cè)線的縱向剖面圖，由于所得結(jié)果是關(guān)于X=0 m對(duì)稱的，因此在這里只展示X=0 m、200 m、400 m、800 m、1 400 m處的剖面圖，視電阻率剖面圖如圖15所示，相位剖面圖如圖16所示。

圖15　視電阻率縱向切片圖Fig．15　The cross-section of apparent resistivity with X=0 m,200 m,400 m,800 m,1 400 m respectively(a)x=0 k m;(b)x=0.2 km;(c)x=0.4 km; (d)x=0.8 km; (e)x=1.4 km

圖16　相位縱向切片圖Fig．16　The cross-section of phase with X=0 m,200 m,400 m,800 m,1 400 m respectively(a)x=0 km;(b)x=0.2 km;(c)x=0.4 km; (d)x=0.8 km; (e)x=1.4 km

取X=0 km處的二維模型進(jìn)行TE模式正演，三維時(shí)X=0 km處的剖面圖與二維TE模式剖面圖對(duì)比,如圖17、圖18所示。

從圖17及圖18可看出，三維X=0 km處視電阻率圖及相位圖都分別與二維TE模式的視電阻率圖和相位圖相似但并不完全相同。需要考慮到三維正演時(shí)用的是三維低阻體，并不是二維模型，而二維本身是對(duì)實(shí)際情況的近似。

4結(jié)論及建議

1)三個(gè)模型的計(jì)算對(duì)比充分說(shuō)明了算法和程序是可靠及有效的。

2)基于MPI實(shí)現(xiàn)了三維大地電磁頻點(diǎn)間的并行，有效地減少了計(jì)算時(shí)間。

3)試驗(yàn)中采用乘大數(shù)法加入第一類邊界條件有時(shí)會(huì)得到錯(cuò)誤結(jié)果，可以進(jìn)一步研究。

4)輔助場(chǎng)的求取對(duì)最后結(jié)果影響較大，如何求得高精度的輔助場(chǎng)很值得研究。

感謝：

譚捍東教授提供的三維模型數(shù)值模擬結(jié)果，在此表示感謝。

圖17　三維X=0 km處剖面圖與二維TE模式視電阻率剖面圖Fig．17　The cross-section of apparent resistivity, the left is the section of 2D with TE model and the right is X=0 km section of 3D(a)二維TE模式視電阻率剖面圖；(b)三維X=0處視電阻率剖面圖

圖18　相位剖面圖Fig．18　The cross-section of phase(a)二維TE模式相位剖面圖；(b)三維X=0處相位剖面圖

參考文獻(xiàn)：

[1]TING S C,HOHMANN G W.Integral equation modeling of three-dimensional magnetotelluric response [J].Geophysics,1981,46(2):182-197.

[2]WANNAMAKER P E.Advances in three-dimensional magnetotelluric modeling using integral equations [J].Geophysics,1991,56(11):1716-1728.

[3]WANNAMAKER P E,HOHMANN G W,WARD S H.Magnetotelluric responses of three-dimensional bodies in layered earths[J].Geophysics,1984,49(9):1517-1533.

[4]MACKIE R L,MADDEN T R,WANNAMAKER P E.Three-dimensional magnetotelluric modeling using difference equations-Theory and comparisons to integral equation solutions [J]. Geophysics,1993,58(2):215-226.

[5]MACKIE R L,SMITH J T,MADDEN T R.Three-dimensional electromagnetic modeling using finite difference equations:The magnetotelluric example[J].Radio Science,1994,29(4):923-935.

[6]譚捍東,余欽范,BOOKER J.大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)報(bào),2003,46(5):705-715.

TAN H D,YU Q F,BOOKER J.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J].Chinese Journal of Geophysics,2003,46(5):705-715.(In Chinese)

[7]MOGI T.Three-dimensional modeling of magnetotelluric data using finite element method[J].Journal of applied geophysics,1996,35(2):185-189.

[8]MITSUHATA Y,UCHIDA T.3D magnetotelluric modeling using the T-Ω Ω finite-element method[J].Geophysics,2004,69(1):108-119.

[9]PUZYREV V,KOLDAN J,DE LA PUENTE J,et al.A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling[J].Geophysical Journal International,2013,193(2):678-693.

[10]SILVA N V D,MORGAN J V,MACGREGOR L,et al.A finite element multifrontal method for 3D CSEM modeling in the frequency domain[J].GEOPHYSICS,2012,77(2):E101-E115.

[11]顏小洋,張偉文,布社輝,等.基于MPI/OPENMP混合編程的三維粒子模擬并行優(yōu)化[J].華南理工大學(xué)學(xué)報(bào)：自然科學(xué)版,2012,40(4):71-78.

YAN X Y,ZHANG W W,BU S H,et al.Parallel Optimization of Three-Dimension Particle Simulation Based on Mixed MPI/OPENMP Programming[J].Journal of South China University of Technology (Natural Science Edition),2012,40(4):71-78.(In Chinese)

[12]童孝忠.大地電磁測(cè)深有限單元法正演與混合遺傳算法正則化反演研究[D].長(zhǎng)沙：中南大學(xué),2008.

TONG X Z.Research of forward using finite method and regularized inversion using hybrid genetic algorithm in Magnetotellurics sounding[D].Changsha:Central South University,2008.(In Chinese)

[13]JIN J.The finite element method in electromagnetics[M].John Wiley & Sons,2014.

收稿日期：2015-04-13改回日期：2015-07-23

基金項(xiàng)目:國(guó)家科技973項(xiàng)目(2014CB440905)；國(guó)家重點(diǎn)實(shí)驗(yàn)室“十二五”項(xiàng)目群(SKLODG-ZY125-01)

作者簡(jiǎn)介：馬駒(1990-)，男，碩士，主要從事大地電磁數(shù)值模擬研究，E-mail:maju_cug08@163.com。 *通信作者：肖調(diào)杰(1991-),男,碩士,從事大地電磁數(shù)值模擬，E-mail:xiaotiaojie@mail.gyig.ac.cn。

文章編號(hào)：1001-1749(2016)03-0289-08

中圖分類號(hào)：P 631.3

文獻(xiàn)標(biāo)志碼：A

DOI：10.3969/j.issn.1001-1749.2016.03.01

Three dimensional modeling of magnetotelluric data using finite element method and MPI parallel computation

MA Ju1，2, XIAO Tiao-jie1，2*, WANG Yun1,3

(1. Institute of Geochemistry, Chinese Academy of Sciences, Guiyang550002, China;2.University of Chinese Academy of Sciences,Beijing100049,China;3. Institute of Geology and Geophysics, Chinese Academy of Science, Beijing100029, China)

Abstract:This paper deals with the forward modeling of three-dimensional (3D) magnetotelluric problem by using the finite element method with fortran language. The validity of this FEM approach is investigated and confirmed by comparing modeling results with a layered medium model, a 2D prism model and a 3D model. The bi-conjugate gradient method is applied in solving the governed equations, which reduced the dimensions of the linear system under the first boundary condition.In addition, the memory requirement is reduced obviously in which only nonzero elements stored. The solution of the linear system is obtained by using a massive parallel program based on MPI between frequencies which massively reduced the computing time. A 3D model has been calculated, and the speed-up ratio is up to 2.328 when 4 processes used.

Key words:magnetotelluric; three-dimensional; modeling; parallel programing; MPI