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

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

?

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

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

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

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

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

0引言

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

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

1三維正演基本理論

1.1邊值問題

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

(1)

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

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

1.2六面體剖分

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

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

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

(2)

其中：ξi、ηi、γi是節(jié)點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)

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

1.4張量阻抗計算

(6)

(7)

(8)

(9)

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

(10)

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

2模型計算

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

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

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

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

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

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

2.2二維棱柱體模型

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

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

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

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

圖8　XY模式與TM模式對比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。計算頻點為3 Hz、1 Hz及0.1 Hz共3個頻點，將YX模式、XY模式分別與參考文獻(xiàn)[6]中的數(shù)值模擬結(jié)果進(jìn)行了對比，YX模式對比如圖10所示，可以看到不管是XY模式還是YX模式，3 Hz及1 Hz時二者擬合很好， 0.1 Hz時曲線形態(tài)基本一致，值有所差別，但最大相對誤差都在5%以內(nèi)。

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

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

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

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

3MPI并行計算

3.1MPI并行實現(xiàn)

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

PROGRAM MT3DForwardFEM

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

……

!初始化MPI系統(tǒng)；得到所有參加運算的進(jìn)程的個數(shù)，放NUMPROCS中；得到當(dāng)前正在運行的進(jìn)程的標(biā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)

……

!為每個進(jìn)程都分配長度一樣的一維數(shù)組用來存儲每個進(jìn)程得到的結(jié)果，F(xiàn)N是頻點數(shù)，!UNX是測點數(shù)，在這里統(tǒng)一分配(FN/NUMPROCS+1)*UNX大小，后面對結(jié)果再進(jìn)行一下處理，收集結(jié)果時就可以用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)用子程序計算形成總體系數(shù)矩陣、加入邊界條件值、解方程、計算視電阻率及相位

… …

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()

……

寫入結(jié)果

……

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

END PROGRAM MT3DFowardFEM

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

3.2并行效率提升評價

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

對圖13、圖14及表1進(jìn)行分析：當(dāng)開啟進(jìn)程數(shù)增加時，并行效率降低；隨著進(jìn)程數(shù)量的增加，并行加速比繼續(xù)增大，但并行效率反而下降，這是因為隨著進(jìn)程數(shù)的增加，用于進(jìn)程間通信的時間所占比例逐漸增大；1個進(jìn)程的并行計算時間比串行時間長，這是因為采用并行需要額外的開銷，總的看來，并行效果還是很明顯的，當(dāng)開啟4個進(jìn)程時達(dá)到了2.3的加速比。

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

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

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

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

由于并行計算并沒有改變算法，因此并行結(jié)果和串行結(jié)果一樣。取X=-1 400 m、-800 m、-400 m、-200 m、0 m、200 m、400 m、800 m、1 400 m 測線的縱向剖面圖，由于所得結(jié)果是關(guān)于X=0 m對稱的，因此在這里只展示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模式正演，三維時X=0 km處的剖面圖與二維TE模式剖面圖對比,如圖17、圖18所示。

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

4結(jié)論及建議

1)三個模型的計算對比充分說明了算法和程序是可靠及有效的。

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

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

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

感謝：

譚捍東教授提供的三維模型數(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.大地電磁法三維交錯采樣有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)報,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é)報：自然科學(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]童孝忠.大地電磁測深有限單元法正演與混合遺傳算法正則化反演研究[D].長沙：中南大學(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

基金項目:國家科技973項目(2014CB440905)；國家重點實驗室“十二五”項目群(SKLODG-ZY125-01)

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

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

中圖分類號：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