陳曉潔，鮑獻豐，李瀚宇，張愛清，周海京

（1.中國工程物理研究院 高性能數(shù)值模擬軟件中心，北京 100088；2.北京應用物理與計算數(shù)學研究所，北京 100094）

飛機、艦船等平臺電磁環(huán)境效應仿真，以及大型陣列天線和天線罩輻射特性仿真等問題均面臨計算模型電大尺寸多尺度瓶頸[1?3]，而大規(guī)模并行技術是解決這一瓶頸問題的最有效手段[4?5].與有限元和矩量法相比，時域有限差分（finite difference time domain，F(xiàn)DTD）方法在時域和空域上均采用顯式迭代，其計算規(guī)模和計算能力不受解法器限制，具備天然的高并行可擴展性，因此大規(guī)模并行FDTD 在實際工程中應用廣泛[6?7].

然而，在實際工程應用中，受到計算資源和成本限制，需要盡可能降低內(nèi)存需求和計算時間，實現(xiàn)“高效”求解.對于電大多尺度模型，采用非均勻網(wǎng)格剖分可以降低內(nèi)存需求.但是，由于穩(wěn)定性條件的限制，F(xiàn)DTD 算法的時間步長取值由最小網(wǎng)格尺寸決定，導致針對多尺度模型的時間步數(shù)較多，計算時間較長.無條件穩(wěn)定FDTD 方法可使得時間步長取值不受穩(wěn)定性條件的限制，但是隱格式無條件穩(wěn)定FDTD方法需要求解大規(guī)模矩陣方程[8]，顯格式無條件穩(wěn)定FDTD 方法需要求解特征矩陣方程[9]，二者均失去了傳統(tǒng)FDTD 方法高并行可擴展性的優(yōu)勢，難以實現(xiàn)大規(guī)模并行計算.另一方面，亞網(wǎng)格FDTD 方法采用粗網(wǎng)格剖分整個計算區(qū)域，并在包含精細結構的區(qū)域采用細化網(wǎng)格進行局部剖分.與基于帶狀非均勻網(wǎng)格剖分的FDTD 方法相比，其不但可以降低內(nèi)存需求，還可以在粗網(wǎng)格區(qū)域和細化網(wǎng)格區(qū)域分別采用大、小時間步長迭代，且不需要求解矩陣方程或特征矩陣方程，在提高計算效率的同時，還保留了傳統(tǒng)FDTD 方法的高并行可擴展性優(yōu)勢，因此并行亞網(wǎng)格FDTD 方法是高效解決電大多尺度模型電磁仿真問題的有效手段.

本文針對電大多尺度模型“高效”電磁仿真需求，基于對稱半正定亞網(wǎng)格FDTD 算法[10]，通過提出融合JASMIN 并行框架[11]h-自適應時間積分算法流程的亞網(wǎng)格FDTD 算法流程，實現(xiàn)了具備大規(guī)模并行能力的亞網(wǎng)格FDTD 算法.該算法采用粗、細兩個網(wǎng)格層分別存儲整個計算區(qū)域網(wǎng)格和局部區(qū)域細化網(wǎng)格，并將兩個網(wǎng)格層均衡地劃分為多個網(wǎng)格片，以網(wǎng)格片為單元實現(xiàn)線程?進程兩級并行.同時，算法支持粗、細網(wǎng)格層分別采用大、小時間步進行時間步進.數(shù)值計算結果表明，并行亞網(wǎng)格FDTD 算法具有較高的精度和并行可擴展性，與并行均勻網(wǎng)格FDTD算法計算結果相比，其計算效率可提高數(shù)10 倍，具備廣闊的實際應用前景.

1 對稱半正定亞網(wǎng)格FDTD 算法

本文簡要介紹對稱半正定亞網(wǎng)格FDTD 算法[10]的原理.分別采用下標c和f定義粗網(wǎng)格和細網(wǎng)格場量，下標i和b定義粗細網(wǎng)格內(nèi)部和粗細網(wǎng)格邊界場量，即eci和ecb分別表示粗網(wǎng)格內(nèi)部電場和粗細網(wǎng)格邊界處粗網(wǎng)格電場，efi和efb分別表示細網(wǎng)格內(nèi)部電場和粗細網(wǎng)格邊界處細網(wǎng)格電場，其位置關系如圖1所示.

圖1 亞網(wǎng)格上場量的位置關系Fig.1 Position relation of field quantities in the sub-grid region

定義電場向量

其中Ne=#eci+#ecb+#ecf，Ne,sub=#eci+#efb+#ecf.定義插值算子矩陣P滿足

將式(1)和(2)代入式(3)，則插值算子P可寫成如下形式

式中，矩陣Pfc反映了粗細網(wǎng)格邊界處粗網(wǎng)格電場和細網(wǎng)格電場的插值關系，即

將時域Maxwell 方程中的電場和磁場分量對空間偏導數(shù)采用中心差分格式離散，電場方程兩邊對時間t求偏導數(shù)后將其代入磁場方程，并基于式(4)對插值算子的定義，可得到如下形式的亞網(wǎng)格FDTD算法矩陣方程[10]：

式中，Dε和Dσ分別為介電常數(shù)和電導率形成的對角矩陣，

式中：Dμc和Dμf分別為粗、細網(wǎng)格磁導率形成的對角矩陣；Sh,c和Se,c與 粗網(wǎng)格區(qū)域空間步長相關；Sh,f和Se,f與細網(wǎng)格區(qū)域空間步長相關.由FDTD 算法的穩(wěn)定性條件可知，若式(6)的矩陣Ctotal為對稱半正定矩陣，則其特征值全部為實數(shù)，選擇適當?shù)臅r間步長，算法將滿足收斂性條件.由于Dε為 對角矩陣，只要Stotal為對稱半正定矩陣，則Ctotal必為對稱半正定矩陣.因此，只要插值算子P的選擇可以保證Stotal為對稱半正定矩陣，則該插值方法必是穩(wěn)定的.

基于上述分析，對稱半正定亞網(wǎng)格FDTD 算法在粗細網(wǎng)格層內(nèi)部均采用傳統(tǒng)的FDTD 方法實現(xiàn)電場和磁場迭代，在粗細網(wǎng)格層邊界需采用以下兩個步驟實現(xiàn)粗、細網(wǎng)格層電場數(shù)據(jù)的插值和交換.

步驟1：粗細網(wǎng)格邊界處粗網(wǎng)格到細網(wǎng)格的電場數(shù)據(jù)細化.

粗細網(wǎng)格邊界處，細網(wǎng)格電場包含兩部分，一部分位于與粗網(wǎng)格電場所在邊重合的位置，記為efb,e，另一部分位于粗網(wǎng)格面上，記為efb,f.如圖2 (a)所示粗細網(wǎng)格交界面處，在ecb已 知的情況下，efb,e可直接取與其所在邊重合位置的ecb，efb,f可基于ecb,1與ecb,2線性插值計算得到.因此，矩陣Pfc第i列為

圖2 粗細網(wǎng)格邊界處電磁場分布示例Fig.2 Example of the electromagnetic field distribution at the boundary of the coarse and fine grids

步驟2：細網(wǎng)格層到粗網(wǎng)格層的電場數(shù)據(jù)粗化.

在粗細網(wǎng)格邊界上，粗網(wǎng)格電場ecb的迭代求解同時涉及粗、細網(wǎng)格上的磁場分量，如圖2 (b)所示.因此，可將ecb分解為兩部分

式中：ecb,p為粗網(wǎng)格磁場的貢獻部分，可基于粗網(wǎng)格磁 場hc1和hc2進 行 迭代求 解；ecb,m為細 網(wǎng)格磁場的 貢獻部分，其求解涉及兩個過程，一是基于細網(wǎng)格磁場插 值得到 圖2(b)所示hf1和hf2，二是基 于hf1和hf2進 行迭代求解.第一個插值過程涉及式(8)右端第二項(PTP)?1PTSf f P，其采用的插值方法是保證算法穩(wěn)定性的關鍵.若直接將步驟1 的矩陣P代入式(8)，難以保證Stotal是對稱半正定矩陣.因此，為了確保穩(wěn)定性，需要基于式(11)，對式(8)的具體實現(xiàn)進行適當修正，使其滿足對稱半正定條件.具體修正過程可參見文獻[10]，在此不再贅述.

2 大規(guī)模并行實現(xiàn)

本文基于北京應用物理與計算數(shù)學研究所自研JASMIN 框架[11]實現(xiàn)亞網(wǎng)格FDTD 算法的大規(guī)模并行計算.JASMIN 框架軟件體系結構如圖3 所示，包括并行計算支撐層、數(shù)值共性層和應用接口層.JASMIN 框架屏蔽了結構網(wǎng)格的大規(guī)模并行和網(wǎng)格自適應計算技術，支持應用軟件開發(fā)人員“并行思考、串行編程”地研制高效使用超級計算機的超級并行應用軟件.

圖3 JASMIN 框架軟件體系結構Fig.3 Software architecture of JASMIN framework

JASMIN 框架的自適應結構網(wǎng)格由嵌套的多個網(wǎng)格層構成，每個網(wǎng)格層被劃分為多個網(wǎng)格片，物理量均存儲于網(wǎng)格片上，以網(wǎng)格片為單元實現(xiàn)線程-進程兩級并行.基于JASMIN 框架實現(xiàn)傳統(tǒng)FDTD 算法時，只需要定義單個網(wǎng)格層.而亞網(wǎng)格FDTD 算法則需要定義兩個網(wǎng)格層，分別為粗網(wǎng)格層和細網(wǎng)格層，前者包含整個計算區(qū)域，后者只包含精細結構區(qū)域.粗網(wǎng)格層和細網(wǎng)格層分別被劃分為多個網(wǎng)格片，這些網(wǎng)格片被分配到各個CPU 進行獨立計算，基于網(wǎng)格片存儲的物理量數(shù)據(jù)稱為數(shù)據(jù)片.

對于只包含單個網(wǎng)格層的傳統(tǒng)FDTD 算法，JASMIN 框架只需負責在網(wǎng)格片之間執(zhí)行數(shù)據(jù)通信操作；而對于亞網(wǎng)格FDTD 算法，JASMIN 框架除了負責在粗、細網(wǎng)格層同層執(zhí)行數(shù)據(jù)通信操作外，還需負責粗、細網(wǎng)格層之間的數(shù)據(jù)傳輸，這一操作流程由h-自適應時間積分算法完成.h-自適應時間積分算法包含時間同步和時間細化兩種模式.前者粗細網(wǎng)格層場量均采用細網(wǎng)格層時間步迭代，后者粗網(wǎng)格層場量采用大時間步迭代，細網(wǎng)格層電磁場采用小時間步迭代.由于時間細化模式的求解效率遠高于時間同步模式，因此本文的計算方法和計算結果均基于時間細化模式.JASMIN 框架屏蔽了h-自適應時間積分算法復雜的數(shù)據(jù)傳輸實現(xiàn)細節(jié)，僅開放了時間、空間細化和空間粗化3 個插值算子接口.理想情況下，亞網(wǎng)格FDTD 算法在粗、細網(wǎng)格層內(nèi)的FDTD 迭代計算流程與經(jīng)典FDTD 算法相同，唯一不同之處在于需基于JASMIN 框架提供的接口實現(xiàn)時間插值算子，以及粗化和細化兩個空間插值算子，以完成粗、細網(wǎng)格層之間的數(shù)據(jù)傳輸.然而，如式（12）所示，為了保證算法的穩(wěn)定性，對稱半正定亞網(wǎng)格FDTD 方法中粗細網(wǎng)格邊界處粗網(wǎng)格電場的求解需同時涉及粗、細網(wǎng)格層的磁場，JASMIN 框架提供的粗化插值接口無法滿足實際計算需求，必須對FDTD 迭代流程進行整體設計.下面將詳細介紹本文基于時間細化模式設計的對稱半正定亞網(wǎng)格FDTD 算法迭代流程.

假設粗細網(wǎng)格細化率為M，粗網(wǎng)格層時間步長為 ?t，則在時間細化模式下，粗網(wǎng)格層電磁場每迭代一步，細網(wǎng)格層電磁場將以時間步長 ?t/M迭代M步.為了論述方便，在粗網(wǎng)格層上定義內(nèi)部電場eci和磁場hc，在細網(wǎng)格層上定義內(nèi)部電場efi和磁場hf，在粗細網(wǎng)格邊界上定義ecb、ecb,p、ecb,m和efb.為了實現(xiàn)時間插值和空間插值，ecb、ecb,p和ecb,m在粗、細網(wǎng)格層均需開辟內(nèi)存.采用ecb（粗）、ecb,p（粗）和ecb,m（粗）表示存儲于粗網(wǎng)格層的場量，采用ecb（細）、ecb,p（細）和ecb,m（細）表示存儲于細網(wǎng)格層的場量.

基于時間細化模式的亞網(wǎng)格FDTD 算法實現(xiàn)流程如圖4 所示.算法1、2 和3 分別具體給出圖4 所示粗網(wǎng)格層電磁場迭代一個時間步、細網(wǎng)格層電磁場迭代M個時間步，以及粗化3 個步驟的具體實現(xiàn)方法.

圖4 亞網(wǎng)格FDTD 算法實現(xiàn)流程Fig.4 Implementation process of the sub-grid FDTD method

算法1 粗網(wǎng)格層電磁場在第n個時間步的迭代流程

1.粗網(wǎng)格層磁場更新

2.粗網(wǎng)格內(nèi)部電場更新

3.在粗細網(wǎng)格邊界上，粗網(wǎng)格的電場更新

算法1 實現(xiàn)了粗網(wǎng)格層電磁場在第n個時間步的迭代流程.標注2～3 如下，算法2 包含M個時間步迭代流程；其中第1 步和第2 步代碼分別基于h-自適應時間積分算法的時間插值接口和空間細化插值接口開發(fā)，以完成粗網(wǎng)格層數(shù)據(jù)到細網(wǎng)格層數(shù)據(jù)的映射；第3 步～第7 步實現(xiàn)了細網(wǎng)格層電磁場的迭代流程.算法3 中第1 步代碼基于h-自適應時間積分算法的空間粗化插值接口開發(fā)，以完成細網(wǎng)格層數(shù)據(jù)到粗網(wǎng)格層數(shù)據(jù)的映射；第2 步實現(xiàn)粗細網(wǎng)格邊界粗網(wǎng)格電場的完整計算.算法2 中第4 步和第7 步實現(xiàn)算法分別為本文第1 節(jié)所示對稱半正定亞網(wǎng)格FDTD 算法實現(xiàn)步驟1 粗細網(wǎng)格邊界粗網(wǎng)格電場到細網(wǎng)格電場的細化和步驟2 細網(wǎng)格層到粗網(wǎng)格層的電場數(shù)據(jù)粗化.

算法2 細網(wǎng)格層電磁場在粗網(wǎng)格層第n個時間步的M步迭代流程

3 數(shù)值算例

3.1 算例1：三維自由空間電磁波傳播問題

自由空間區(qū)域的左下角和右上角直角坐標分別為(?1.8, ?1.8, ?1.8) m 和(4.3, 4.3, 4.3) m，其x、y和z方向的空間網(wǎng)格步長均為0.1 m.在空間區(qū)域 [(1, 1,1) m, (1.4, 1.4, 1.4) m] 內(nèi)采用亞網(wǎng)格剖分，網(wǎng)格細化率為20，即亞網(wǎng)格區(qū)域x、y和z方向的空間網(wǎng)格步長均為0.005 m.入射平面波沿x軸傳播，極化方向沿z軸，其電場分量表達式為Einc=z?2(t?t0?x/c)e(t?t0?x/c)2/τ2，c為光速，取值3×108m/s，τ=2×10?8s，t0=4τ.觀察點P1和P2的坐標分別為(0.1, 1.3, 1.3)和 (1.3, 1.3, 1.3)，分別位于亞網(wǎng)格區(qū)域外和亞網(wǎng)格區(qū)域內(nèi).

圖5 為采用本文的亞網(wǎng)格FDTD 算法的4 核并行計算結果與理論結果的比對，可以看出二者吻合較好.

圖5 監(jiān)測點時域電場波形Fig.5 Time domain electric field waveform of the monitoring points

3.2 算例2：多模導彈電磁耦合問題

多模導彈CAD 模型如圖6 所示.該模型彈體長度約為6 m，導引頭包含天線罩、拋物面天線和線纜等精細結構，天線罩厚度為1.5 mm，其中天線罩為聚四氟乙烯材料，其他部件為理想導體.入射平面波沿-z軸傳播，極化方向沿x軸，時域信號源為調(diào)制高斯脈沖，其中心頻率為1.5 GHz，帶寬為3 GHz.包含多模導彈模型的FDTD 計算區(qū)域為[(?0.505, ?0.505,?0.505) m，(0.5, 0.5, 5.6) m].采用亞網(wǎng)格FDTD 方法計算時，粗網(wǎng)格層網(wǎng)格步長為7.5 mm，網(wǎng)格索引范圍為[(0, 0, 0), (133, 133, 813)]，網(wǎng)格總數(shù)約為14.6×106；細網(wǎng)格層覆蓋整個導引頭區(qū)域，其計算區(qū)域?qū)诖志W(wǎng)格層的網(wǎng)格索引范圍為[(42, 42, 688), (92, 92, 747)]，網(wǎng)格步長為1.5 mm，即網(wǎng)格細化率為5，細網(wǎng)格層總網(wǎng)格數(shù)約為19.5×106，因此粗細網(wǎng)格層總網(wǎng)格數(shù)約為34.1×106.導引頭區(qū)域的粗、細網(wǎng)格層計算網(wǎng)格模型如圖7 所示.由于導引頭區(qū)域包含的天線罩、線纜等精細結構，粗網(wǎng)格層無法實現(xiàn)幾何模型的精確剖分，出現(xiàn)了結構破損、缺失等現(xiàn)象，引入局部加密的細網(wǎng)格層后則有效避免了上述現(xiàn)象.

圖6 多模導彈CAD 模型Fig.6 The CAD model of the multi-mode missile

圖7 導引頭區(qū)域的網(wǎng)格模型Fig.7 The grid model of the seeker area

為了驗證亞網(wǎng)格FDTD 算法的正確性和有效性，采用全局加密的均勻網(wǎng)格FDTD 算法與亞網(wǎng)格FDTD算法的計算結果進行比對.均勻網(wǎng)格FDTD 算法的空間步長與亞網(wǎng)格FDTD 算法細網(wǎng)格層的空間步長一致，為1.5 mm，其總網(wǎng)格數(shù)約為998.9×106.兩種方法總計算物理時間均為29 ns，亞網(wǎng)格FDTD 算法中，粗網(wǎng)格層計算時間步數(shù)約為2 000 步，細網(wǎng)格層時間步數(shù)約為10 000 步，均勻網(wǎng)格FDTD 算法的計算時間步為10 000 步.

圖8 給出了亞網(wǎng)格FDTD 方法及全局加密均勻網(wǎng)格FDTD 方法計算得到的相同物理時刻導引頭附近時域磁場分布偽彩圖，可見二者吻合較好；圖9 給出了導引頭內(nèi)部監(jiān)測點的時域電場波形，兩種方法計算得到的電場波形吻合較好，由此驗證了亞網(wǎng)格FDTD 方法的計算精度和高效性.

圖8 多模導彈導引頭附近時域場分布（物理時刻4.3 ns）Fig.8 Time domain field distribution near the missile seeker of the multi-mode (physical time is 4.3 ns)

圖9 多模導彈導引頭內(nèi)部監(jiān)測點時域電場波形Fig.9 Time domain electric field waveform of the internal monitoring points of the multi-mode missile seeker

表1 給出了兩種方法的計算開銷和加速比，其中加速比Sp=T1/Tp，T1為基于1 個CPU 核的計算耗時，Tp為基于p個CPU 核的計算耗時.在相同并行規(guī)模下（128 CPU 核并行），亞網(wǎng)格方法內(nèi)存占用僅為全局加密方法的7%，計算速度提高54.5 倍.由此可見，對于多尺度模型，亞網(wǎng)格FDTD 方法具有較大的性能優(yōu)勢.

表1 亞網(wǎng)格FDTD 方法與均勻網(wǎng)格FDTD 方法計算開銷對比Tab.1 Comparison of the computational overhead between the sub-grid FDTD method and the uniform grid FDTD method

為驗證亞網(wǎng)格FDTD 方法的并行可擴展性，基于多模彈模型開展了強可擴展并行效率測試（即計算規(guī)模保持不變的前提下，改變并行規(guī)模），測試得到的并行效率和加速比數(shù)據(jù)如表2 所示.對于強并行可擴展性測試，隨著并行規(guī)模的增加，通信耗時占比越來越大.表2 中，雖然當CPU 核數(shù)從8 增長到16 時，通信耗時已變成主要矛盾，計算效率陡降，但是128 CPU 核相對于串行計算其強可擴展并行效率仍可達到35%，驗證了算法的高并行可擴展性.

表2 亞網(wǎng)格FDTD 方法并行效率和加速比測試結果Tab.2 Testing results of the parallel efficiency and acceleration ratio of the subgrid FDTD method

4 結 論

本文基于高性能并行計算框架JASMIN，實現(xiàn)了具備大規(guī)模并行能力的對稱半正定亞網(wǎng)格FDTD 算法.在保證精度和穩(wěn)定性的前提下，該算法可支持粗網(wǎng)格區(qū)域和細網(wǎng)格區(qū)域分別采用大時間步長和小時間步長進行迭代計算，且具備高并行可擴展性，可高效地解決電大多尺度模型的電磁模擬問題.數(shù)值算例表明，針對電大多尺度模型，本文提出的并行亞網(wǎng)格FDTD 方法可大大減少內(nèi)存需求，提高計算效率，且具有較高的并行可擴展性，因此具有廣闊的應用前景.