基于GPU的FDTD麥克斯韋方程快速求解_邵楨

1、第 36 卷第 24 期計 算 機 工 程2010 年 12 月Vol .36No. 24Computer EngineeringDecember 2010開發(fā)研究與設(shè)計技術(shù)文章編號：10003428(2010)24027803文獻標(biāo)識碼：A中圖分類號：TP391.9基于 GPU 的 FDTD 麥克斯韋方程快速求解邵楨 a，蔡紅星 b，徐春鳳 a(長春理工大學(xué) a. 計算機科學(xué)與技術(shù)學(xué)院；b. 理學(xué)院，長春 130022)摘 要：采用圖形處理器(GPU)為主計算核心，應(yīng)用時域有限差分法(FDTD)實現(xiàn)電磁學(xué)中麥克斯韋方程組的快速求解。通過對 FDTD 求解麥克斯韋旋度方程的直接時間域的分析，給

2、出 FDTD 的仿真算法。根據(jù) GPU 能高效地提高 FDTD 的仿真速度，解決 FDTD 仿真算法中的計算量龐大問題。利用 GPU 在 FDTD 計算中的處理能力，實現(xiàn)了更長的脈沖持續(xù)時間和龐大的模型求解與仿真，在適當(dāng)?shù)臅r間內(nèi)完成了超大量的仿真計算。根據(jù)在 CPU 和 FDTD 上的實際計算結(jié)果表明，基于 GPU 的 FDTD 仿真算法具有高精度和高效率等特點。關(guān)鍵詞：時域有限差分法；圖形處理器；麥克斯韋方程；電磁場仿真Solving Quickly of Finite Difference Time Domain Maxwells Equations Based on GPUSHAO Zh

3、ena, CAI Hong-xingb, XU Chun-fenga(a. School of Computer Science and Technology; b. School of Science, Changchun University of Science and Technology,Changchun 130022, China)【Abstract】Using Graphics Processing Units(GPU) as main computational core, the Finite Difference Time Domain(FDTD) is presente

4、d for solving electromagnetic Maxwells equations. The FDTD simulation algorithm is developed by analyzing the direct time-domain solution like FDTD Maxwells curl equations. It is analyzed and compared with GPU that how GPU can be used to greatly speed up FDTD simulations. So enormous computation pro

5、blem in FDTD simulations is resolved. Using GPU processing power in FDTD calculations, it implements much longer pulse lengths and larger models, and completes expensive simulations computation in reasonable time. According to the comparing of actual computation results in CPUs and FDTD, It is prove

6、d that FDTD simulations based on GPU is accurate and efficient.【Key words】Finite Difference Time Domain(FDTD); Graphics Processing Units(GPU); Maxwells equations; electromagnetic field simulation1 概述時域有限差分法(Finite Difference Time Domain(FDTD) 是求解電磁學(xué)中麥克斯韋方程式的重要方法之一，一直以來獲得了廣泛的使用1。它能夠求解任何介質(zhì)材料組成的任意幾何圖形的

7、模型。另外，F(xiàn)DTD 還能很容易地處理較短的電磁波的脈沖信號( 函數(shù))1，使其有利于解決電磁學(xué)中的許多問題。但是，對于像 FDTD 這樣大模型空間的仿真或者較長的非正弦波信號的處理，都需要超大量的浮點計算，甚至在巨型高性能計算系統(tǒng)中都要運行幾個月或者更多的時間。即使是具有單元數(shù)量級約為 109 的小空間模型，對于一個毫秒級持續(xù)時間波形的仿真，傳統(tǒng)的高性能計算系統(tǒng)的 1 000 個處理器也要計算至少一個月的時間2。所以，基于傳統(tǒng) CPU 編寫的現(xiàn)階段的 FDTD 程序?qū)τ诮鉀Q該類問題是非常困難的，需要進一步研究提高這些計算仿真速度的方法。2 用于通用計算的 GPU圖形處理器(Graphics P

8、rocessing Units, GPU)具有低成本、高性能的特點，是一種向量計算進行了優(yōu)化的高度并行的數(shù)據(jù)流處理機2。顯性結(jié)構(gòu)的數(shù)據(jù)并行使 GPU 不再需要大量復(fù)雜的控制邏輯，數(shù)據(jù)流經(jīng)高速內(nèi)存接口時不再需要大量的緩存。這樣就可以空出大量的空間用于存儲計算單元。最近幾年，GPU 已經(jīng)被廣泛地使用到絕大多數(shù)的通用計算中。目前，GPU 在解決 FDTD 這樣需要大規(guī)模計算量的問題上顯現(xiàn)出了巨大的潛能。278在浮點運算方面，GPU 的良好性能大大高于 CPU，歸功于 GPU 具有專用的流架構(gòu)。近年來，GPU 的可編程性和浮點運算的高精度，使其能夠出色地完成絕大多數(shù)的計算工作。在原始數(shù)據(jù)浮點計算方面，

9、GPU 的性能就遠高于傳統(tǒng)的通用CPUs，并且具有更高的內(nèi)存帶寬。如 GeForce 8800 GTX 的浮點運算性能達到峰值為 345 GFLOPS，內(nèi)存接口的速度為 86.4 Gb/s，而英特爾酷睿 2 主頻為 2.66 GHz(1 GHz 的前端總線 ) 四核 CPU 處理器的浮點運算的峰值理論上僅為 21.3 GFLOPS 和 10.7 Gb/s 的內(nèi)存帶寬2；另外 GPU 性能提升的速度也是遠快于 CPU，GPU 在 6 個月的時間里性能可以提升 1 倍，而 CPUs 需要 18 個月。GeForce 8800 GTX 顯卡具有 128 個流處理器，數(shù)據(jù)的并行處理只需要計算流而不再需

10、要共享內(nèi)存空間，分支邏輯在硬件低開銷的情況下具有了更高的效能。因為采用流處理的數(shù)據(jù)不再需要駐留內(nèi)存，這樣每個指令需要大量計算，其算法就顯得尤為重要。而許多傳統(tǒng)的算法是不能直接運行在流架構(gòu)上的，只有采用并行處理、具有大量計算、高速的算法才可以在流架構(gòu)上運行。基于現(xiàn)有的 GPU 產(chǎn)品，本文編寫了新的 FDTD 仿真程序，可以處理麥克斯韋方程的求解和大型模型空間的仿真，其速作者簡介：邵 楨(1974)，女，講師、碩士，主研方向：計算機仿真，智能算法；蔡紅星，副教授、博士；徐春鳳，講師收稿日期：2010-09-16E-mail：這就是最初由 K.S.Yee 提 出 的 時 域

11、有限差分法 FDTD，值得一提的是其著名的 Yee 空間網(wǎng)格模型，電場 E 和磁場 H 各節(jié)點的空間排布如圖 1 所示，又稱為 Yee 氏元胞圖4。度得到了較大的提升。3 時域有限差分(FDTD)FDTD 用于求解任意模型空間的麥克斯韋方程，此外， FDTD 還可以幫助人們求解那些相對難度很大或者用分析法不太可能實現(xiàn)的模型。FDTD 求解麥克斯韋旋度方程的直接時間域的方法如下：仿真的時間分辨力，波形和信號源的持續(xù)時間決定了時間步數(shù)。問題空間以三維網(wǎng)格的形式存儲，網(wǎng)格的每個單元設(shè)計成有相應(yīng)介電性能的材料類型，存儲初始化為 0 的電場 E 和磁場 H 的 x、y、z 方向的分量。初始化的空間采用的

12、基本 FDTD 算法如圖 2 所示。 E = H(1)t H = E+(2)Jt首先把麥克斯韋旋度方程式 (1)和式(2) 中的失量轉(zhuǎn)換成想 x、y、z 3 個方向的分量，然后通過中心差分法近似地將任意分量表示成時間和空間的一階微分3。電場強度矢量 E 的 x 分量如式(3)所示，其他 5 個分量的表示方法相似：Ex=1HzH ytyzEy=1 HxHztzxEz=1H yHxtxyHx=1EyEztzyH y=1EzExtxzHz=1ExEytyxExEyEz(3)(4)(5)(6)(7)圖 2FDTD 的頂層算法在該算法中，大部分時間消耗在向量 E 和 H 修正計算的(8) 多重循環(huán)上，每

13、個網(wǎng)格空間的 E 向量修正和 H 向量修正計算是獨立的并且可以并行計算的。算法執(zhí)行的復(fù)雜度和內(nèi)存需求估算如下：圖 1三維 Yee 氏元胞圖由圖 1 可見，每一個磁場分量由 4 個電場分量環(huán)繞；同樣，每一個電場分量由 4 個磁場分量環(huán)繞。電場和磁場在時間順序上交替抽樣，抽樣時間彼此相差半個時間步，使麥克斯韋旋度方程離散以后構(gòu)成顯式差分方程，從而可以在時間上進行迭代求解。給定相應(yīng)電磁問題的初始值和邊界條件以后，利用 FDTD 遞進公式就可以遞進地求得各個時刻各單元節(jié)點處的場值。4FDTD 基本算法FDTD 算法在解麥克斯韋方程時首先是以時間步 n 來執(zhí)行每個網(wǎng)格的 E 向量修正方程，然后以時間步

14、n+1/2 來執(zhí)行每個網(wǎng)格的 H 向量修正方程5。模型的空間解析能力決定了算法復(fù)雜度 時間步數(shù) 網(wǎng)格數(shù) 6 操作數(shù)其中，因數(shù) 6 代表在 3D 模型中需要修正的分量的個數(shù)。分量的修正計算正如前面所述，可以簡化為非均質(zhì)三維模型空間的每個網(wǎng)格分量的 21 位浮點操作。內(nèi)存需求 網(wǎng)格數(shù) 6 常數(shù)其中，因數(shù) 6 代表分量數(shù)；常量表示一個浮點數(shù)所占的字節(jié)數(shù)，1 個字節(jié)常用來表示網(wǎng)格的材料類型。5 在 GPU 上的實現(xiàn)與數(shù)據(jù)分析FDTD 修正計算既有數(shù)據(jù)并行又有高強度計算，這樣就非常適合在 GPU 上執(zhí)行。滿足數(shù)據(jù)并行條件的模型空間的每個網(wǎng)格都要執(zhí)行相同的向量 E 和 H 的修正計算，每個修正計算至少每秒

15、執(zhí)行 18 個浮點操作，這樣足以忽略大多數(shù) GPU 的內(nèi)存延遲問題。最初的 FDTD 系統(tǒng)由 C 和 Cg 編寫，應(yīng)用于 GeForce 7 系列的具有特殊的頂點和碎片處理器的 GPU 上。存儲在存儲器中的向量 E 和 H 作為 32 位浮點三原色的二維紋理。三原色 RGB 通道存放 x、y、z 3 個方向的分量。50%的精密亮度紋理以指針流的形式存儲模型空間的材料類型。在這 2 種情況下，三維空間被夷平為二維紋理，通過三維到二維的地址轉(zhuǎn)換來訪問空間的某一點。這就使得整個三維空間可以在一個渲染遍中更新，避免了在寫入數(shù)據(jù)損壞后的非法讀取。材料類型的指針流可以用于查找存儲在紋理中的介質(zhì)材料屬性。E

16、 和 H 散射場的更新計算轉(zhuǎn)換到由 Cg 編寫的碎片程序來進行。著色器更新程序把存儲在紋理中向量 E 和 H 作為輸入數(shù)據(jù)，把數(shù)據(jù)渲染輸出到幀緩沖對象 FBO 中。每個著色器更新程序的計算是以渲染一個四邊形作為初始化設(shè)置的。279利用 Linux、C#、OpenGL 和 GeForce 7 系列 GPU，本文編寫了一套適用于流處理器的 FDTD 法求解麥克斯韋方程的仿真程序 。 圖形硬件可以通過標(biāo)準(zhǔn)的 OpenGL 來存取 。 OpenGL 紋理先將 FDTD 模型空間轉(zhuǎn)換到 GPU 設(shè)備自身的內(nèi)存中，再通過 OpenGL 2.0 應(yīng)用程序接口(API)將其經(jīng)由幀緩沖器傳送到主機中。在 FDT

17、D 更新計算中產(chǎn)生的碎片可以利用 GPU 的碎片處理器來處理。在測試中，穩(wěn)定極限為 0.9，一個 020 GHz 的點源置于網(wǎng)格中心，如圖 3 所示。仿真程序運行環(huán)境分別是 Intel P4 2.4 GHz、512 MB 內(nèi)存和 256 MB 顯存的 NVIDIAs GeForce 7800 GS。圖 7 256256 網(wǎng)格性能曲線圖 8 512512 網(wǎng)格性能曲線圖 3 仿真范圍求解圖 4 是一個時間步函數(shù)的高速仿真精度曲線圖。修正計算的過程是 8 個8 192 個時間步。平均累計誤差的速度是CPU 的 12.943 9 倍。圖 4 精確度對比圖 5圖 10 是一個時間步函數(shù)在不同網(wǎng)格數(shù)情況

18、下，GPU和 CPU 執(zhí)行仿真程序的性能曲線圖。由圖中數(shù)據(jù)可計算出隨著仿真規(guī)模的擴大，GPU 下執(zhí)行的仿真程序平均速度是以 CPU 速度的 17.533 7 倍的比率進行遞增。圖 5 6464 網(wǎng)格性能曲線圖 6 128128 網(wǎng)格性能曲線280圖 9 10241024 網(wǎng)格性能曲線654 CPU32 GPU102 0484 0968 1921 024時間步圖 10 20482048 網(wǎng)格性能曲線6 結(jié)束語對于 FDTD 這樣的具有高強度計算量并且計算數(shù)據(jù)流只有有限的或沒有數(shù)據(jù)相關(guān)性的高性能的計算系統(tǒng)，GPU 的使用具有較好的前途。用 C 和 Cg 編寫的 GPU FDTD 程序表現(xiàn)出了比運行

19、在 CPU 上的傳統(tǒng) FDTD 程序加速許多倍的性能。在通常情況下，如果模型空間非常大以至于導(dǎo)致碎片處理器總是處于忙碌狀態(tài)，那么 GPU 程序也比 CPU 程序快幾倍，即使是最慢的七系列的 GPU，如 GeForce Go 7400。如果模型空間不是足夠大，不足以使碎片處理器飽和，那么 GPU 程序的執(zhí)行速度反要比 CPU 程序慢。因此，如果 FDTD 模型空間足夠大，就可以更好地發(fā)揮 GPU 效能。當(dāng)模型單元數(shù)的數(shù)量級是 1403 以上時，GPU 處理 FDTD 更新計算的速度是傳統(tǒng)多核 CPU(例如 Intel 或者 AMD)處理速度的 12 倍之多。另一個有趣的現(xiàn)象是在 GPU 上進行少

20、量的迭代計算時初始化流的過程是有效能損失的，但是隨著迭代次數(shù)的增加，這個初始化損失的效能就被均攤了。因為典型的 FDTD 仿真算法需要成千上萬次的迭代才能計算出結(jié)果，所以除了非常小的模型空間，在大多數(shù)情況下，使用 GPU 的速度要遠遠快于使用 CPU 產(chǎn)品的速度。本文應(yīng)用 GPU 平臺和 FDTD 方法，較好地解決了以麥克斯韋方程組為代表的具有邊界條件復(fù)雜、計算數(shù)據(jù)量大、數(shù)據(jù)相關(guān)性低等特點的方程組的快速求解問題，該工作將促進此類方程求解技術(shù)的發(fā)展。(下轉(zhuǎn)第 283 頁)則返回(8)。(10) 在點 ( xi , yi , zi ) 上添加顏色值 cx ，如果 color ( cx ). CRC

21、apability Course ( si ).StudentNum ， 并 且 c x Cii ， 則DCT ( xi , yi , z i ) = cx ， Course ( si ).CountInWeek 1 。 cx 添加到顏色集合 Cii 中。(11)返回(5)。(12)遍歷 z 軸的課程編號，如果存在任意課程編號 z j ， Course ( z j ).CountInWeek 0 ( z j )，則返回(6)；否則排課完畢。3 結(jié)果與分析根據(jù)上述排課表模型(每周 5 天，每天 5 個課次)和優(yōu)化模型，在教師、班級、課程確定的情況下，可以給出教室數(shù)量 的 范 圍 NumberRa

22、nge 。 設(shè) 定 NumberRange 的 上 限 是 MaxNum ，下限是 MinNum ，并且教師數(shù)量為 Tn ，班級數(shù)量為 Cn ，一周內(nèi)總課程次數(shù)為 Sn ，第 i 個教室一周使用次數(shù)為 CTi 。由于現(xiàn)實中教師數(shù)量總是大于班級數(shù)量，因此，約定 Tn Cn 。(1) 首先給出 MinNum 的值 。 當(dāng) NumberRange = MinNum時，對于任意 i ，1 i MinNum ，CTi =25?？傻茫篗inNum Sn25 = Sn ，則 MinNum =。(2)然后給出 MaxNum 的值。假定每個教師只教授 1 個班級，每個班級每次課用一個教室。當(dāng) NumberRan

23、ge = MaxNum 時，因為 Tn Cn ，所以此時 Cn 個班級同時上課，則 MaxNum = Cn 。實際上，每個教師可能教授多個班級，設(shè)定 m 個教師教授 1 個班級， n 個教師教授 2 個班級， n 個教師教授的班級編號為 C11 , C12 , C1n , C 21 , C 22 , C2n 。 對于 C11 與C21 , C12 與 C22 , C1n 與 C2n 共 n 個組合，每個組合的 2 個班級有且只有一個上課，因此 MaxNum = Cn n ，得到：Sn NumberRange Cn - n25用計算機實現(xiàn)本校電信學(xué)院 2007 年-2008 年第 1 學(xué)年排課表

24、，共有 64 個任課教師，45 個班級，37 門課程，一周總課次為 196，使用 8 個教室(最少教室數(shù)目)，得到表 1。在表 1 中，每個單元格有若干二元組。二元組的第 1 個數(shù)字表示課程編號，一個課程編號表示該課程的一個課次，第 2 個數(shù)字表示該課次所在的教室編號。結(jié)果滿足基本約束條件、硬約束條件及部分軟約束條件。表 1 結(jié)果是在教室數(shù)最少的情況下得到，如果教室數(shù)是既定的，則可得到最優(yōu)排課表。文獻4中應(yīng)用“關(guān)聯(lián)規(guī)則”解決排課表的沖突問題，文獻5通過“矩陣解法”解決排課表問題，它們都是針對排課表本身。本文給出的模型主要目的之一是給出教學(xué)資源消耗的最小值，并盡可能使問題的解決模式化，這是與傳統(tǒng)排課表模型的最大區(qū)別。表 1排課表結(jié)果周一周二周三周四周五第 1 節(jié)第 2 節(jié)第 3節(jié)第 4節(jié)第 5節(jié)4 結(jié)束語排課表問題是一個經(jīng)典問題，本文在分析和總結(jié)以往排課表研究成果的基礎(chǔ)上，提出基于立方體染色的排課表模型。在該模型上引入立方體染色的定義，設(shè)計排課表優(yōu)化模型，很好地解決了排課表問題。傳統(tǒng)排課表模型通常只關(guān)注排課問題的最優(yōu)解，而本文最大的區(qū)別是將雜亂的排課表的約束條件轉(zhuǎn)化為立方體的點線面的關(guān)系。參考文獻1 廖飛雄, 馬 良. 圖著色問題的啟發(fā)式搜索螞蟻算法J. 計算機工程, 2007, 33(16): 191-192, 195.2 Bo