艦船磁特征磁矩量法的圖形處理單元加速計算研究

郭成豹,劉大明

（海軍工程大學電氣工程學院,湖北武漢 430033）

艦船磁特征磁矩量法的圖形處理單元加速計算研究

郭成豹,劉大明

（海軍工程大學電氣工程學院,湖北武漢 430033）

磁矩量法（MMM）非常適合用于艦船等鐵磁目標的磁場特征建模,因為這種方法不需要對空氣等非鐵磁區(qū)域進行網格劃分,然而這種方法需要存儲和處理稠密的滿陣,計算時間很長,矩陣存儲需要非常大的存儲空間。聯合采用MMM和多層自適應交叉近似（MLACA）法雖然可以顯著縮減內存需求和計算時間,但對于求解精細劃分網格的問題,計算時間還是太長。將具有強大并行計算能力的圖形處理單元（GPU）用于加速具有天然良好并行特性的MLACA法,并給出了相應的并行計算格式,實現了艦船磁特征MMM的大規(guī)模并行加速計算。典型算例結果表明,GPU并行計算的加速比超過120倍,對于精細劃分為100 000薄殼單元的艦船殼體,其計算時間也僅有約4.3 min.采用該方法的計算結果與商業(yè)有限元軟件相比差別小于1%,為艦船磁場的大規(guī)模建模提供了一種快速、精確、簡便的數值計算工具。

電磁學；艦船磁場；磁矩量法；多層自適應交叉近似；圖形處理單元

0 引言

對艦船、汽車等三維復雜鐵磁目標的磁場特征進行快速準確分析,已經廣泛應用于艦船磁場探測、艦船消磁、交通控制等諸多領域[1-4]。郭成豹等曾經聯合采用磁矩量法（MMM）和多層自適應交叉近似（MLACA）法進行了艦船磁場特征的精確建模分析[3],但是計算時間非常長（5萬個單元計算時間約6.5 h）,在實際使用中非常不方便。

MLACA法具有天然的并行化特性,可以采用并行化技術進行加速。圖形處理單元（GPU）由于具備強大的并行處理能力和極高的存儲器帶寬,可以用于科學運算、數據分析、線性代數等需要大量重復數據運算和密集內存存取的場合。計算統(tǒng)一設備架構（CUDA）是美國NVIDIA公司于2007年推出的GPU通用計算產品,是第一種直接使用類C語言進行通用計算的開發(fā)環(huán)境和軟件體系,不需借助圖形學應用程序接口。由于在性能、成本和開發(fā)時間上比傳統(tǒng)的中央處理單元（CPU）解決方案具有顯著優(yōu)勢,CUDA已經在地球物理、流體力學、信號探測等領域獲得廣泛應用,并取得了豐碩的成果[5]。

本文將GPU并行計算技術引入艦船磁場特征建模計算中,所采用的GPU硬件為NVIDIA C2075,擁有448個CUDA核心,6GB GDDR5顯存。將艦船劃分為一定數目的三角形薄殼單元,分別采用GPU和CPU進行了感應磁場特征計算。計算結果表明, GPU的計算結果與CPU計算結果一致,但極大地縮短了計算時間。采用本文技術編制了高性能計算軟件,可以使得上述5萬單元的模型計算時間縮短為約2 min,加速比超過120倍。即使對于精細劃分為10萬單元的復雜結構艦船,其計算時間也僅有約4.3 min,非常便于實用。

1 磁矩量法

艦船鐵磁殼體處于磁化場Hc（地磁場）中,由于磁化場幅值較小,薄板的相對磁導率μr可以認為是線性的。由于相對磁導率μr（對于造船鋼板,一般μr＞100）,薄板內的磁化強度可以認為是沿著薄板的切向分布。艦船殼體可以劃分為N個單元,每個單元內部的磁化強度假定是均勻的。采用配點法,可以得到磁介質內以M為變量的積分方程

式中：HcSi是在鋼板面上的磁化磁場；e是薄板厚度； Lj是源單元j的邊界；rij是場單元i到源單元j的矢量；rij是場單元i到源單元j的距離。將磁介質區(qū)域劃分為n個單元,當單元足夠小時,就可以認為單元內的磁化強度M和磁化率χ均為恒定值,于是（1）式即可轉化為方程組

2 多層自適應交叉近似法

一個鐵磁目標可被劃分為很多單元,這些單元又可按照結構關系歸類為若干單元組。根據上述單元組之間的相互關系,MMM的整個系數矩陣可以劃分為一系列的子矩陣。對角線上的子矩陣描述了距離相隔較近的單元組之間的相互作用,可直接采用標準MMM計算。非對角線上的子矩陣是非滿秩的,描述的是距離相隔較遠的單元組之間的相互作用。

距離較遠的源單元組Nn與場單元組Nm之間的子矩陣Zmn是非滿秩的,可以采用ACA法壓縮為幾個行向量uk和列向量vk的乘積[6-8]式中：k是矩陣中向量序號。因此,相對于全填充的系數矩陣,這種壓縮技術可以節(jié)省內存,加快計算速度。為了停止近似過程,本文采用了下述停止準則,

式中：‖‖F表示Frobenius范數；Sk是當前低秩近似式,其秩＜k；γ是停止門限值。對于良好分離的單元組,壓縮子矩陣的秩遠遠小于源單元組變量數目和場單元組變量數目。

為了獲得更高的矩陣壓縮比,需采用MLACA法。按照二叉樹的模式將鐵磁目標空間遞歸劃分為小的區(qū)域,使得各個小區(qū)域內的單元組符合如下要求：

式中：B（N）是包含N的最小球半徑；dist（Nn,Nm）是單元組中心之間的距離；η≥1是確定單元組良好分離與否的參數。MLACA法的內存需求趨近于O（nlogn）,計算費用趨近于O（n2）.一艘典型艦艇采用MLACA法計算得到的壓縮系數矩陣如圖1所示,圖中數字表示壓縮子矩陣的秩。

3 算法集成和GPU并行化

為了使得軟件集成化,采用了模塊化的集成框架[9],如圖2所示。具體計算過程如下。

3.1 多層區(qū)域分解

將艦船劃分為若干三角形薄殼單元,根據各單元重心點的坐標將單元進行多層區(qū)域分解。此部分計算耗時非常短,10萬單元只要幾秒鐘的時間,因此可以不用并行化。

3.2 MMM系數積分和MLACA法計算及其GPU加速

圖1 壓縮矩陣的圖形表示Fig.1 Graphical representation of the compressed matrix

圖2 MMM和MLACA法的集成框圖Fig.2 Block diagram of the integration framework with MMM and MLACA

分析MLACA算法的計算流程可以知道,主要的計算時間是系數矩陣行或列的計算。因此,可以將矩陣行和列的計算在GPU中并行進行,計算完后返回CPU,每個矩陣單元一個線程。ACA算法中的很多矢量點乘、范數、絕對值、搜索最大值等操作,可以在CPU上執(zhí)行,這是因為此類操作數據量巨大而操作較為簡單,在GPU上進行并無優(yōu)勢。

系數積分計算會得到4個子矩陣塊：Zxx、Zxy、Zyx、Zyy,考慮到MLACA算法是一種純代數方法,為了更有效地降低計算和存儲復雜度,將4個子矩陣塊作為一個整體矩陣Z進行壓縮處理。同時,每次進行更多系數矩陣元素的計算可以減小GPU的數據通信開銷。

系數積分的GPU并行計算過程如下。

輸入：

開始時輸入到GPU顯存,以后不再重復輸入,直接調用即可,可以提高數據訪問速度。

輸出：

MMM的積分系數子矩陣Zmn.

步驟1并行計算所有單元的重心坐標{gex,

步驟2并行計算所有單元的坐標轉換矩陣

步驟3GPU的線程th1,…,thn分別計算Zm1, Zm2,…,Zmn,得到Zmn={Zm1,Zm2,…,Zmn}.每次計算系數矩陣的行或列時,從主機端向GPU傳輸的僅是單元序號集合,從而減少數據傳輸時間。

3.3 GMRES法求解及其GPU加速

由于易于并行化且收斂速度快,Saad等在1986年提出的廣義最小殘差法（GMRES）非常適用于線性方程組ZM=Hc的求解[10]。該方法通過Arnoldi迭代產生Krylov子空間K（Z,q1,k）的一組標準正交基K=span{q1,Zq1,…,Zk-1q1},并以此正交基表示{Mk},使Mk=M0+span{q1,q2,…,qk},并使得‖Hc-ZMk‖達到極小,Mk即為方程組的近似解。GMRES法迭代計算過程中,90%以上的計算量是線性代數運算,運算時間主要消耗在具有較高復雜度的矩陣向量乘法運算上。因此,本文將CUDA用于矩陣向量乘法的并行化運算。采用GPU作為協(xié)處理器實現系統(tǒng)方程的并行求解,CPU則負責算法的迭代流程控制。為了提高運算效率,在GPU上將矩陣和向量分組,用一個線程執(zhí)行對應的一組乘法操作,多個線程并行執(zhí)行。由于得到的矩陣性態(tài)較為良好,本文的求解計算并未采用預條件技術。

4 計算實例及結果

本文計算的硬件環(huán)境如下：CPU為Intel Xeon CPU 2687 W,內存為64GB DDR3- 1600；GPU為NVIDIA Tesla C2075,顯存為6GB GDDR5存儲器。軟件環(huán)境為Microsoft Visual Studio 2010,NVIDIA CUDA 5.5.算例的計算分析流程如圖3所示,采用三角形薄殼單元,MLACA法和GMRES法的收斂判據都設為10-3.

算例1和算例2為直徑10 m、厚度0.05 m、相對磁導率160的空心球殼,分別劃分為2 048個和8 192個三角形薄殼單元；算例3為長32 m、寬5 m、厚度0.05 m、相對磁導率160的潛艇殼體,劃分為99 098個三角形薄殼單元。3個算例均處于外界均勻磁化磁場Bz=50 000 nT作用下。

圖3 算例的計算分析流程Fig.3 Computation flowchart of the examples

算例的各部分耗時、加速比等如表1所示,從表中可以看出,基于GPU的CUDA并行計算相對于CPU串行計算,對于稀疏剖分（2 048個單元）的空心球殼算例1和精細剖分（8 192個單元）的空心球殼算例2,總加速比均超過120倍。對于算例3,潛艇殼體劃分為接近10萬單元,GPU并行計算的求解時間約為4.3 min,速度非?？?而相應的CPU計算時間預計將達數十小時,在此沒有進行計算實驗。此外,算例的GPU并行計算結果與CPU串行計算結果也進行了對比,二者差別小于0.5%,可忽略不計。

算例2和算例3計算得到的空心球殼和潛艇殼體的磁化強度分布分別如圖4和圖5所示,圖示結果與商用有限元軟件計算結果具有非常高的吻合度,差別小于1%.

為了驗證本文所提出方法的磁場計算精度,在潛艇殼體中心線正下方設置了一條磁場計算線（如圖6所示）,深度10 m,長度50 m,分別采用本文提出的方法和商業(yè)有限元軟件計算上述計算線上的潛艇殼體磁場,如圖7所示,發(fā)現二者計算結果差別亦小于1%.因此,本文提出的方法對于潛艇殼體外部磁場的計算精度非常高。

表1 算例初始條件、耗時和加速比Tab.1 The initial conditions,computation times,speedup ratios of the examples

圖4 球殼的磁化強度計算值Fig.4 The calculated values of magnetization of spherical shell

圖5 潛艇殼體的磁化強度計算值Fig.5 The calculated values of magnetization of submarine hull

圖6 潛艇殼體磁場計算線Fig.6 The evaluating line of magnetic fields from the submarine hull

圖7 潛艇殼體磁場的商業(yè)有限元軟件計算結果與本文方法計算結果的對比Fig.7 The comparison of the magnetic fields from the submarine hull calculated by the FEM commercial software and the proposed method

本文的系數矩陣計算采用雙精度,但是計算結果采用單精度存儲,因為單精度數據存儲量僅有雙精度的50%,可以大幅度節(jié)省存儲空間。為了進一步減小數值誤差,本文計算中Z矩陣系數后一致乘以磁化率χ,χ=μr-1,以改善積分方程矩陣的性態(tài),并在相應的求解結果里乘以磁化率χ求得真實結果。

5 結論

本文采用GPU并行計算技術對基于MLACA法和MMM的艦船磁場特征建模過程進行了加速。典型算例表明獲得了超過120倍的加速比,并且計算結果與商業(yè)有限元軟件相比差別小于1%,從而使得艦船磁場的MMM大規(guī)模建模技術實用化。即使精細劃分為10萬個三角形薄殼單元的艦船殼體,其計算時間也僅有約4.3 min,非常便于實用。隨著GPU性能的不斷改善和可編程性的不斷提高,GPU法將在艦船磁場特征計算分析中發(fā)揮更大的作用。

（References）

[1] 郭成豹,肖昌漢,劉大明.基于積分方程法和奇異值分解的磁性目標磁場延拓技術研究[J].物理學報,2008,57（7）： 4182-4188.

GUO Cheng-bao,XIAO Chang-han,LIU Da-ming.Research on the continuations of magnetic field of magnetic object based on integral equation method and singular value decomposition[J].Acta Physica Sinica,2008,57（7）：4182-4188.（in Chinese）

[2] Chadebec O,Coulomb J,Janet F.A review of magnetostatic moment method[J].IEEE Transactions on Magnetics,2006,42（4）： 515-520.

[3] 郭成豹,劉大明.薄鋼殼物體磁特征建模研究[J].兵工學報,2012,33（8）：912-915.

GUO Cheng-bao,LIU Da-ming.Modeling of magnetic signatures of thin sheet objects[J].Acta Armamentarii,2012,33（8）：912-915.（in Chinese）

[4] Nguyen T S,Guichon J M,Chadebec O,et al.Ships magnetic anomaly computation with integral equation and fast multipole method[J]. IEEE Transactions on Magnetics,2011,47（5）：1414-1417.

[5] Sanders J,Kandrot E.CUDA by example：an introduction to general-purpose GPU programming[M].Boston：Addison-Wesley Professional,2010.

[6] Zhao K,Vouvakis M N,Lee J F.The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems[J].IEEE Transactions on Electromagnetic Compatibility,2005,47（4）：763-773.

[7] Bebendorf M.Approximation of boundary element matrices[J]. Numerische Mathematik,2000,86（4）：565-589.

[8] Katrijn F,Marc V B.Solving a large dense linear system by adaptive cross approximation[J].Journal of Computational and Applied Mathematics,2010,234（11）：3181-3195.

[9] Hugo G E,Yaniv B,Amir B,et al.Software framework for integration of method of moments kernels with the adaptive cross approximation（ACA）algorithm[C]∥2010 IEEE 26th Convention of Electrical and Electronics Engineers in Israel.Eilat,Israel： IEEE,2010：634-636.

[10] Saad Y,Schultz M.GMRES：A generalized minimal residual algorithm for solving nonsymmetric linear systems[J].SIAM J Sci Stat Comput,1986,7（3）：856-869.

GPU-accelerated Magnetic Moment Method for Magnetic Signatures of Ship

GUO Cheng-bao,LIU Da-ming
（School of Electrical and Information Engineering,Naval University of Engineering,Wuhan 430033,Hubei,China）

Magnetic moment methods（MMMs）are particularly well suited to model the magnetic signatures of the ferromagnetic objects,such as ship,because they do not require the meshing of non-ferromagnetic materials like air.However the methods need the storing and solving of a fully dense matrix system.Computation takes a long time,and a large storage space is reqiured for storing the matrix.The method based on MMM and multilevel adaptive cross approximation（MLACA）leads to a significant reduction of the required memory and computation times.But the computation times are still too long for the resolution of very fine meshed geometries.The graphical processing units（GPU）with high parallel processing capability is used to accelerate MLACA with the natural parallelism,and the parallel computation strategies are presented,implementing the fast analysis of large scale MMM for magnetic signatures of ships.The results of typical examples indicate that the speedup ratio of the parallel computation of graph-ics processing unit（GPU）is more than 120,and the computation time for the ship hull finely meshed to 100 000 sheet elements is only about 4.3 min.The difference between the calculated results of the proposed method and the FEM commercial software is less than 1%.Thus,it provides a fast,accurate and simple method for the large scale modeling of ship magnetic signatures.

electromagnetism；ship magnetic field；magnetic moment method；multilevel adaptive cross approximation；graphics processing unit

U665.18

A

1000-1093（2014）10-1638-06

10.3969/j.issn.1000-1093.2014.10.018

2014-01-03

國家自然科學基金項目（51277176）

郭成豹（1975—）,男,講師,博士。E-mail：guochengbao@outlook.com；劉大明（1956—）,男,教授,博士生導師。E-mail：liudaming@126.com