并行時域有限差分法網格自動剖分技術

摘要： 為了構建適合于時域有限差分法求解的離散幾何模型，擺脫大量繁雜的手工操作，提出了一種高效的大規模并行時域有限差分法網格自動剖分技術.該并行剖分方法將三角面元計算機輔助設計模型文件作為輸入數據，能夠自動建立包含多種介質的任意三維實體網格模型，且其并行執行過程中無需過多的數據交換，并行效率可達99%.基于該離散模型自動生成技術，采用時域有限差分方法，計算了含微帶板的簡易計算機機箱耦合效應，其計算結果與商用電磁仿真軟件的計算結果相吻合，表明該剖分方法的準確性.最后，通過250個處理器核，并行剖分了網格規模約6億的某真實計算機機箱模型，其并行剖分時間僅為0.2 s，驗證了該并行剖分方法的高效性，表明該網格剖分技術能有效地解決大規模并行時域有限差分法的離散幾何建模問題.

關鍵詞： 時域有限差分法；網格剖分；并行計算；計算機輔助設計

中圖分類號： O441.4文獻標志碼： ACartesian Mesh Generator for Parallel

FiniteDifference TimeDomain MethodCHEN Linglu1，ZHOU Haijing2，LI Hanyu2，FU Haijun1，3，LIAO Cheng1

（1. Institute of Electromagnetics， Southwest Jiaotong University， Chengdu 610031， China； 2. Institute of Applied Physics and Computational Mathematics， Beijing 100094， China； 3. Beijing Institute of Electronic System Engineering， Beijing 100854， China）

Abstract：In order to generate discrete geometricmodel that can be solved by the finitedifference timedomain （FDTD） method， and avoid complicated manual operations， an efficient parallel mesh generation technique for massively parallel FDTD simulation was presented. The parallel mesh generator is capable of dissecting arbitrary 3D multimaterial models， with trianglefacet computer aided design files as input data. Without the need for excessive data communication between the processors， its parallel efficiency was close to 100% in the test. The coupling effect of a simple computer case with a microstrip board was simulated by the FDTD method， based on this mesh generator. The results by the FDTD coincided with the ones calculated by the commercial electromagnetic simulation software， demonstrating the correctness of this mesh dissection method. Finally， when dissecting a real computer case with 250 processors and the grid number of 0.6 billion， the CPU time was only about 0.2 s. This efficient mesh dissection indicates that the proposed parallel mesh generation technique can efficiently solve the modeling in massively parallel FDTD simulation.

Key words：finitedifference timedomain （FDTD）； mesh generation； parallelization； computer aided design （CAD）

時域有限差分法（finite difference time domain， FDTD）是一種簡單、靈活的電磁仿真算法，且其易于并行，現已被各界廣泛應用與發展[15].為了使FDTD仿真算法走向實用化，高效的自動幾何建模方法是必不可少的，這在大規模并行FDTD方法的運用中尤為必要.近年來，越來越多的學者基于計算機輔助設計（computer aided design， CAD）研究FDTD離散幾何建模方法.如文獻[6]對AutoCAD軟件進行二次開發，實現三維實體的FDTD網格自動剖分；本文作者又對這種方法提出了補充和改進[7].但這種方法的剖分過程是在AutoCAD軟件環境下進行的，因此其剖分效果受到AutoCAD軟件的限制，且不易并行化.Y. Srisukh等則提出了一種可以將三角面元數據轉換為可用于FDTD方法計算的六面體網格模型的射線追蹤方法[8]；之后Jeff T. MacGillivray利用該方法在一個單核工作站上完成了網格規模為1萬億的FDTD電磁目標建模，但其剖分時間長達10 h[9]；此外，與之類似的還有文獻[10]提出的投影求交方法等.雖然上述方法都具有較好的魯棒性，但需要花費較多的時間進行求交以及排序運算.為了進一步提高FDTD自動幾何建模的效率，本文將原鑄造業所采用的切片線掃描法[11]引入到電磁幾何建模中，可最小化求交和排序運算.然后，基于并行自適應結構網格應用支撐軟件框架（J parallel adaptive structured mesh applications infrastructure， JASMIN）實現了剖分算法的并行化，以解決單臺計算機內存受限的問題，并加快網格剖分速度.該方法適合于處理結構復雜，網格規模大的大規模并行FDTD仿真計算問題.西南交通大學學報第48卷第4期陳伶璐等：并行時域有限差分法網格自動剖分技術1基于FDTD的網格自動剖分算法及其并行化實現1.1三角面元CAD數據文件及介質識別方法三角面元CAD數據文件所存儲的模型信息類似于有限元網格劃分，以三角面片來擬合逼近實體模型的表面輪廓特征.由于兩個不同三角形的3個頂點坐標不可能完全相同，因此，只需給出各個三角形的頂點坐標，模型的形狀也就可以完全確定了.STL （stereolithography）格式文件是一種標準的三角面元CAD數據文件[12]，幾乎所有的CAD系統都提供了STL格式文件的輸入、輸出接口，如Solid Work、Pro/E、Auto CAD等，且當下一些主流的電磁仿真軟件也都提供了它的輸入、輸出接口，如CST、HFSS等.根據數據儲存形式的不同，STL格式文件可以分為二進制碼和ASCII碼兩種形式.相較于二進制碼存儲形式，ASCII碼形式具有更強的可讀性.如圖1所示，采用ASCII碼存儲形式的STL格式文件以“solid component”、“facet normal”、“outer loop”等關鍵詞來區分各種數據的意義.各個實體的信息以“solid”為開頭，以“endsolid”為結束；單個實體段中，每個三角形的信息以“facet”為開頭，以“endfacet”為結束；法向矢量以“normal”為引導；頂點信息以“outerloop”為開頭，以“endloop”為結束；每個頂點的坐標以“vertex”為引導，每個頂點有x、y、z 3個坐標值.本文所研究的網格自動剖分算法以這種ASCII碼形式的STL格式文件作為輸入數據，具有很好的通用性.

solid component1： facet normal outer loopvertex vertex vertex endloopendfacet…endsolid component1：

圖1STL格式文件的ASCII碼存儲形式

Fig.1An example of an ASCII STL file

由圖1可知，STL格式文件只保存有實體的名稱、每個三角形的法向矢量和頂點坐標這3類信息，并不包含實體模型的介質信息.因此，本文通過預先設定的命名規則對各實體進行命名，以實體名來識別不同實體的介質信息.例如，為“002_solid1”時，“002”為介質編號（如：“002”代表銅），“solid1”為自定義的實體名稱.本文所有STL格式文件由商用軟件CST導出，在CST幾何建模平臺中，即可對各實體按上述命名規則進行命名.1.2FDTD網格自動剖分算法本文研究的網格自動剖分算法先以切片的形式得到實體各切割平面的封閉輪廓，再利用線掃描的方法判斷空間中各網格的具體屬性.圖2展示了一中空立方體的網格剖分過程，其中圖2（b）為該中空立方體的三角面元形式.在剖分前，預先設X、Y、Z方向的剖分間距分別為Δx、 Δy 和 Δz，根據FDTD算法的要求，剖分間距應小于λ/12[13]，其中λ為FDTD仿真的最小波長.具體剖分步驟如下：

步驟1選取某個坐標平面（如XOY面）為基準面，再以平行于基準面的切割平面z=（k+1/2）Δz對實體進行切割，其中k為代表第k個切割平面的整數.描述實體表面輪廓的三角面元將與切割平面相交而產生交線.這樣，在切割平面上就可得到一個由若干線段首尾相接而組成的封閉環，如圖2（c）所示.

步驟2在切割平面上以直線y=（j+1/2）Δy（其中，j為整數，代表第j條掃描線）對該平面進行線掃描.這些掃描線將可能與封閉環相交，如圖2（d）所示.

步驟3將掃描線與封閉環的交點按順序編號為1，2，…，n，并按“奇數偶數”配對，如奇數交點1和偶數交點2為一組，奇數交點3和偶數交點4為一組，以此類推，如圖2（e）中粗實線所示，則可以判斷奇數交點與偶數交點之間的區域為實體內部區域.

步驟4以直線x=（i+1/2）Δx（其中，i為整數，代表第i條掃描線）掃描上述“奇偶”配對的直線，如圖2（f）所示，其交點皆位于FDTD的六面體網格中心，判定這些網格的屬性為實體.

（a） 中空立方體模型（b） 中空立方體的三角面元形式示例：z=3.5（c） 切割平面 z=（k+1/2）Δz（d） 掃描線y=（j+1/2）Δy（e） 交點編號和“奇偶”配對（f） 掃描線 x=（i+1/2）Δx圖2網格剖分過程示意

Fig.2Mesh dissection process步驟5重復上述步驟，直到全計算空間被切割掃描完畢.

當剖分目標為多實體組成的多介質模型時，則需要對每個實體重復上述步驟.

在網格剖分過程中，需要注意對奇異點的處理.奇異點是指切割平面或掃描線與三角面元的特殊交點.這些奇異點大致可分為2類：第1類是在用切割平面對實體進行切割時，切割平面與三角面元位于同一平面或只過三角面元其中1個頂點而產生的，這種情況將導致切割平面與某三角面元有3個或1個交點，而非1段交線；第2類是在對切割平面進行線掃描時，掃描線只過封閉環某頂點而產生的，這種情況將導致掃描線與封閉環的交點為奇數，不利于程序對網格屬性的判斷.為使程序順利執行，當出現上述情形時，本文將切割平面或掃描線略微偏移一個小量，以避開奇異點的產生.此外，須使切割平面的偏移方向與三角面元外法線方向相反，以確保切割平面偏移后仍與實體相交.1.3網格自動剖分算法的并行化實現本文所提出的FDTD網格自動剖分算法的并行化是在JASMIN框架下實現的.JASMIN框架是一個基于消息傳遞接口庫 （message passing interface， MPI）的并行框架，可支持數萬至數十萬處理器核的大規模并行計算[14].

基于JASMIN框架下“塊網格片層次結構網格層網格片”的并行區域分解策略，將全計算空間分解為若干個子計算空間，令每個網格片只負責處理其中一個子空間.各子計算空間在獲得STL模型文件數據后，在該子空間內獨立地執行網格自動剖分程序，并只對處于該子空間內的實體目標進行網格剖分.由于本文所提出的網格自動剖分算法在并行剖分過程中并不需要進行MPI通信，即無需數據交互，而是直接將剖分后的離散網格模型用于并行FDTD仿真，因此，這種網格自動剖分算法的并行效率是近乎100%的.

除網格自動剖分算法其本身的高效性外，如何使各進程快速獲得STL模型文件數據也是影響網格剖分效率的關鍵因素.最簡單易行的獲取方式為“各處理器并行讀取”，即所有處理器同時讀取目標文件.然而，隨著待剖分模型復雜程度的增加，STL模型文件的數據量也將增大，所用處理器核的數量也將增多.這將導致在程序運行的初始階段，系統的I/O節點負載迅速增大，且在I/O的負載平衡方面也將難以控制，即，若采用上萬處理器核并行讀取一個25 MB大小的STL模型文件，僅讀取文件的時間將長達十幾分鐘.因此，這種“各處理器并行讀取”的方式并不適用于所有情況.另一種獲取方式為僅由主進程讀取目標文件，然后經關鍵詞識別，保留有用信息，再廣播給其他進程.本文所用的STL模型文件最大約25 MB，而其有用信息僅為KB量級，經關鍵詞識別后廣播有用信息可很大程度地減少系統的I/O節點負載，降低程序對硬件設備的要求.采用這種方式獲取一個25 MB大小的STL模型文件數據，時間為秒級.因此，本文采用由主進程讀取，再廣播給其他進程的方式來實現STL模型文件的快速讀取.2數值算例計算機機箱的電磁耦合效應研究一直是電磁兼容、干擾和防護等領域的重要課題[15]，而計算機機箱復雜結構的精確幾何建模也是對其進行仿真分析的難點之一.本文將基于并行FDTD網格自動剖分技術對含微帶板的簡易計算機機箱模型和某真實計算機機箱進行建模仿真，用以驗證所提出的并行網格自動剖分算法的準確性和高效性.2.1簡易計算機機箱模型網格剖分簡易計算機機箱模型如圖3所示，在主板位置放置一個兩端開路的微帶板，該簡易計算機模型尺寸為170 mm×420 mm×430 mm.為設置激勵源和邊界條件，所取計算空間將大于機箱模型尺寸，統一取260 mm×520 mm×520 mm.在一個4核4 GB 2.13 GHz Intel Xeon Linux工作站上計算該算例.

圖3帶微帶板的簡易計算機機箱模型

Fig.3Simple computer case with microstrip board

表1給出了不同網格規模和處理器核數量的簡易計算機機箱模型剖分時間.需要指出的是，當空間步長取值為Δx=Δy=Δz=0.1 mm 時，所對應的網格規模約為700億（2 600×5 200×5 200），已經超出了該工作站4 GB的內存限制，因此，這里的網格剖分數據并未被保存，只作觀察并行效率用.并行效率[16]的定義為

E（q）=Ts/（Tqq），（1）

式中：Ts為串行執行時間；q為使用的處理器核數；Tq為使用q個處理器核并行執行的時間.

表1簡易計算機機箱網格剖分CPU運行時間

Tab.1CPU time of the simple computer case mesh dissection

計算空間網格規模處理器核CPU運行時間/s260 mm×520 mm×520 mm130×260×26010.701 6092 600×5 200×5 20011 161.1292581.3064291.331

從表1可以看出，采用并行網格自動剖分算法可以大大縮減網格剖分所占用的CPU運行時間.對于網格量約700億規模的計算機機箱模型，采用2個處理器核和4個處理器核剖分的并行效率分別為99.87%和99.64%.

由于該算例計算機機箱結構相對簡單，當空間步長取值為Δx=Δy=Δz=2 mm時，即可滿足FDTD程序[17]的仿真精度要求，此時網格規模約900萬（130×260×260）.激勵信號為一高斯脈沖，設脈寬為0.2 ns，時延為0.16 ns，幅值為1 kV/m.入射波為正面垂直入射的垂直極化平面波.以計算機正面細縫內側中心（130 mm，312 mm，47 mm）、計算機機箱中心（130 mm，260 mm，260 mm）和主板中心（62 mm，262 mm，366 mm）為觀察點，對比FDTD程序和CST微波工作室時域求解器的仿真結果.如圖4～6所示，兩者的計算結果幾乎完全吻合，由此驗證了該網格自動剖分算法能準確地將三角面元CAD模型文件數據轉換為FDTD程序可用的六面體網格模型.

圖4簡易計算機機箱正面細縫中點

電場強度時域圖

Fig.4Electric field in the middle of the silt located

on the front of simple computer case圖5簡易計算機機箱中心點電場強度時域圖

Fig.5Electric field at the center of

simple computer case圖6主板中心點電場強度時域圖

Fig.6Electric field at the center of motherboard

2.2真實計算機機箱網格剖分某真實計算機機箱模型如圖7所示，機箱物理尺寸為250 mm×500 mm×550 mm，箱體表面遍布各種形狀的散熱孔，機箱背部有三芯電源線，機箱內有復位線、主板（以介質敷銅板替代）.在曙光5 000高性能計算機上通過并行FDTD程序計算該真實計算機機箱的電磁耦合效應.空間步長取值為Δx=Δy=Δz=0.5 mm，所對應的網格規模約為6億（520×1 020×1 120）.

（a） 實物圖（b） FDTD網格模型圖7某真實計算機機箱

Fig.7A real computer case表2給出了該真實計算機機箱模型的網格剖分時間，在250個處理器核并行剖分的條件下，CPU運行時間僅為0.2 s.FDTD程序的激勵設置與2.1節相同.

表2某真實計算機機箱網格剖分CPU運行時間

Tab.2CPU time for the mesh dissection of the real computer case

計算空間網格規模處理器核/個CPU運行時間/s260 mm×510 mm×560 mm520×1 020×1 1202500.203 658

圖8為該計算機機箱中心點處的時域電場強度曲線.該算例表明了本文提出的并行FDTD網格自動剖分方法的高效性.

圖8某真實計算機機箱中心點電場強度時域圖

Fig.8Electric field at the center point of

the real computer case3結論本文提出了一種基于大規模并行FDTD的網格自動剖分方法，可以將包含多種介質的任意三維實體快速、準確地離散為適用于FDTD計算的六面體網格結構.由于并行剖分過程不需要進行MPI通信，剖分算法的并行效率可達99%.簡易計算機機箱模型的FDTD仿真結果和CST軟件計算結果的對比，驗證了該剖分方法的準確性.以250個處理器核耗時僅0.2 s完成約6億網格規模的某真實計算機機箱網格剖分，表明了該并行剖分方法的高效性.這種以三角面元CAD數據文件作為輸入數據的網格自動剖分算法對于大規模并行FDTD的離散化網格建模是十分有效的.參考文獻：[1]YEE K S. Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media[J]. IEEE Transactions on Antennas and Propagation， 1966， 14（3）： 302307.

[2]WATANABE Y， IGARASHI H. Accelerated FDTD analysis of antennas loaded by electric circuits[J]. IEEE Transactions on Antennas and Propagation， 2012， 60（2）： 958963.

[3]楊丹，方劍，廖成，等. PML吸收邊界條件在孔縫耦合模擬中的應用[J]. 西南交通大學學報，2005，40（1）： 108112.

YANG Dan， FANG Jian， LIAO Cheng， et al. Application of PML absorbing boundary condition in numerical simulation of pulses coupling into slot[J]. Journal of Southwest Jiaotong University， 2005， 40（1）： 108112.

[4]鐘選明，廖成. 基于時間反演和單站天線的導體目標微波成像[J]. 西南交通大學學報，2011，46（3）： 451455.

ZHONG Xuanming， LIAO Cheng. Microwave imaging of perfectly conducting cylinder based on time reversal and single antenna[J]. Journal of Southwest Jiaotong University， 2011， 46（3）： 451455.

[5]廖成，衛濤，陳偉. 整數微分進化策略及其在微波成像中的應用[J]. 西南交通大學學報，2007，42（6）： 647652.

LIAO Cheng， WEI Tao， CHENG Wei. Integer coded differential evolution strategy and application to microwave imaging[J]. Journal of Southwest Jiaotong University， 2007， 42（6）： 647652.

[6]張秋菊，王秉中. 時域有限差分法電磁仿真的網格自動剖分[J]. 電子科技大學學報，2007，36（1）： 6669.

ZHANG Qiuju， WANG Bingzhong. Auto dissection of entity with grid for FDTD[J]. Electromagnetic Simulation Journal of University of Electronic Science and Technology of China， 2007， 36（1）： 6669.

[7]CHEN Linglu. Auto dissection of entity with threedimensional network based on FDTD[C]∥2010 AsiaPacific Symposium on Electromagnetic Compatibility. Beijing： IEEE Press， 2010： 920923.

[8]SRISUKH Y， NEHRBASS J， TEIXEIRA F L， et al. An approach for automatic grid generation in threedimensional FDTD simulations of complex geometries[J]. IEEE Transactions on Antennas and Propagation， 2002， 44（4）： 7580.

[9]MACGILLIVRAY J T. Trillion cell CADbased cartesian mesh generator for the finitedifference timedomain method on a singleprocessor 4GB workstation[J]. IEEE Transactions on Antennas and Propagation， 2008， 56（8）： 21872190.

[10]楊麗霞，葛德彪，白劍，等. 三角面元數據模型FDTD網格生成技術[J]. 西安電子科技大學學報，2007，34（2）： 208302.

YANG Lixia， GE Debiao， BAI Jian， et al. A novel FDTD modeling technique based on triangle meshunits of an object[J]. Journal of Xidian University， 2007， 34（2）： 208302.

[11]梁英業，戴挺，趙建新，等. 基于STL的切片線掃描法網格剖分技術[J]. 鑄造，2005，54（10）： 10021005.

LIANG Yingye， DAI Ting， ZHAO Jianxin， et al. Mesh generation using a slice line scanning method based on STL files[J]. Foundry， 2005， 54（10）： 10021005.

[12]ITO Y， NAKAHASHI K. Surface triangulation for polygonal models based on CAD data[J]. International Journal for Numerical Methods in Fluids， 2002， 39（1）： 7596.

[13]葛德彪，閆玉波. 電磁波時域有限差分方法[M]. 2版. 西安：西安電子科技大學出版社，2005： 3839.

[14]莫則堯，張愛清. 并行自適應結構網格應用支撐軟件框架[R]. 北京：北京應用物理與計算數學研究所，2009： 317.

[15]廖成，楊丹，方劍. 不同極化方向電磁脈沖作用于計算機機箱的效應分析[J]. 微波學報，2007，23（4）： 2831.

LIAO Cheng， YANG Dan， FANG Jian. Numerical analysis for EMP having different polarization direction coupling to a computer case[J]. Journal of Microwaves， 2007， 23（4）： 2831.

[16]喻煒，顏宏，金之雁. 并行效率的初步研究[J]. 應用氣象學報，1996，7（1）： 6168.

YU Wei， YAN Hong， JIN Zhiyan. The preliminary study on the efficiency of parallel model computation[J]. Quarterly Journal of Applied Meteorology， 1996， 7（1）： 6168.

[17]CHANG Lei. Parallel FDTD analysis of largescale ridged TEM horn antenna array[C]∥2009 International Conference on Microwave Technology and Computational Electromagnetics. Chengdu： IEEE Press， 2009： 344347.