摘要：時(shí)域有限差分法（FDTD）求解電磁學(xué)中麥克斯韋方程組是科學(xué)與工程計(jì)算中一個(gè)非常重要的算法。通過對(duì)FDTD 求解麥克斯韋旋度方程的直接時(shí)間域的分析，給出其基于多個(gè)GPU 組成異構(gòu)機(jī)群系統(tǒng)上的并行加速算法，用OpenCL、CUDA 和MPI 編程模型實(shí)現(xiàn)了并行程序。在目前的主流NVIDIA 和ATI 的GPU 平臺(tái)上，加速的并行FDTD 程序相對(duì)CPU 串行程序和8 個(gè)CPU 核的MPI 并行程序，分別獲得了超過8 倍和1.5 倍的加速，并在多個(gè)GPU 卡上獲得了接近線性加速的擴(kuò)展性能。

　　引言

　　Maxwell 方程組用數(shù)學(xué)模型揭示了自然界一切宏觀電磁現(xiàn)象所遵循的普遍規(guī)律，一百多年來，人們依據(jù)這組方程將電磁波的研究滲透到各個(gè)領(lǐng)域，應(yīng)用十分廣泛，例如微波、天線、電磁成像、電磁防護(hù)、無線電傳播、導(dǎo)航、雷達(dá)技術(shù)、地下電磁探測、電磁兼容等等。研究電磁場問題歸根結(jié)底便是在特定的邊界條件下解出Maxwell 方程組的解或近似解。隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，已經(jīng)出現(xiàn)解決Maxwell 方程組的多種有效的數(shù)值計(jì)算方法，而近40 年發(fā)展為迅速的電磁場數(shù)值計(jì)算方法是時(shí)域有限差分（FDTD）方法。FDTD方法是1966 年由K. S. Yee 首次提出，經(jīng)過幾十年的發(fā)展，已經(jīng)應(yīng)用于輻射天線的分析、微波器件和導(dǎo)行波結(jié)構(gòu)的研究、散射和雷達(dá)截面計(jì)算、電磁兼容分析、微帶線與PBG 結(jié)構(gòu)仿真等多個(gè)電磁領(lǐng)域中。而在電磁工程中，存在著一些消耗計(jì)算資源巨大的問題，例如大型天線陣列的分析、復(fù)雜結(jié)構(gòu)電磁特性仿真等等，諸如這些大規(guī)模復(fù)雜電磁工程問題，需要的計(jì)算量都異常龐大，包含了海量的浮點(diǎn)計(jì)算，即使在巨型高性能計(jì)算系統(tǒng)中都要運(yùn)行幾個(gè)月或者更多的時(shí)間。即使是具有單元數(shù)量級(jí)約為109 的小空間模型，對(duì)于一個(gè)毫秒級(jí)持續(xù)時(shí)間波形的仿真，傳統(tǒng)的高性能計(jì)算系統(tǒng)的1000 個(gè)處理器也要計(jì)算至少一個(gè)月的時(shí)間。基于傳統(tǒng)CPU 編寫的現(xiàn)階段的FDTD 程序?qū)τ诮鉀Q該類問題是非常困難的，需要進(jìn)一步研究提高這些計(jì)算仿真速度的方法。

　　GPU（Graphic Processing Unit）是個(gè)人計(jì)算機(jī)中的重要部件，其主要功能是加速圖形處理任務(wù)的吞吐率。自從2000年以來，GPU 的浮點(diǎn)運(yùn)算能力以每12 個(gè)月翻一番的速度成指數(shù)式增長。到2008 年底，主流GPU 的浮點(diǎn)峰值和顯存帶寬就已經(jīng)大大超過了同時(shí)代的X86 CPU,相應(yīng)的，單位功耗可提供的浮點(diǎn)峰值也高出CPU 接近一個(gè)數(shù)量級(jí)。近年來，GPU 在雙浮點(diǎn)計(jì)算部件、訪存層次等方面不斷完善，性能、功耗、成本方面的綜合優(yōu)勢越來越明顯，在高性能計(jì)算領(lǐng)域的應(yīng)用前景越來越廣，是實(shí)現(xiàn)Peta/Exascale計(jì)算技術(shù)的重要途徑。目前已有大量的科學(xué)工程因GPGPU 而獲益，相關(guān)的研究工作表明利用GPU 的計(jì)算能力對(duì)求解FDTD 能夠獲得一定的加速效果。本文根據(jù)FDTD 算法和GPGPU 體系結(jié)構(gòu)的特點(diǎn)，提出了一種基于CPU-GPU 異構(gòu)機(jī)群的FDTD 并行計(jì)算方法，研究了內(nèi)存訪問和數(shù)據(jù)通信優(yōu)化策略，在主流GPGPU 平臺(tái)上獲得了較好的計(jì)算加速效果，并且能夠支持多個(gè)GPGPU 卡的近似線性擴(kuò)展。

　　1 FDTD 算法概述

　　FDTD 算法的基本思想是：使用差分網(wǎng)格，將連續(xù)的電磁場問題變?yōu)殡x散系統(tǒng)問題，用各離散點(diǎn)上的數(shù)值解來逼近連續(xù)場域內(nèi)的真實(shí)解。計(jì)算域空間節(jié)點(diǎn)采用Yee 元胞（圖1），每個(gè)磁場分量由四個(gè)電場分量環(huán)繞，同樣每一個(gè)電場分量由四個(gè)磁場分量環(huán)繞，采用蛙跳格式，電場和磁場在時(shí)間順序上采用交錯(cuò)抽樣，這樣使得麥克斯韋旋度方程離散后構(gòu)成顯式差分方程，利用前一時(shí)刻已知的電場和磁場求出當(dāng)前時(shí)刻的電磁場，在時(shí)間上迭代求解。FDTD采用吸收邊界條件的方法，把整個(gè)計(jì)算域劃分成包括散射體的總場區(qū)以及只有反射波的散射場區(qū)，使得計(jì)算可以在有限的空間范圍內(nèi)進(jìn)行。因此，只要給定的相應(yīng)電磁問題的初始值，通過FDTD 就可以迭代求出此后各個(gè)時(shí)刻空間電磁場的分布，并在時(shí)間軸上逐步推進(jìn)地求解電磁場。

　　利用中心差分法，可以將 MaxWell 標(biāo)量方程中更新某個(gè)時(shí)刻的E 沿Y 軸方向場值寫成下面的顯示遞推方程：

　　另外 5 個(gè)遞推方程也可以類似推導(dǎo)出來，每個(gè)網(wǎng)格點(diǎn)上的場分量的新值依賴于該點(diǎn)在前一時(shí)間步長時(shí)刻的值及該點(diǎn)周圍的臨近點(diǎn)上另一場量在早半個(gè)時(shí)間步長時(shí)的值。

　　這樣，當(dāng)某一時(shí)刻所有場域內(nèi)的電磁場計(jì)算完后，使用迭代過程就可以計(jì)算出隨時(shí)間變化的網(wǎng)格內(nèi)電磁場分布情況。FDTD 算法流程如圖2 所示。

　　為了確定 FDTD 程序中需要加速的函數(shù)，我們使用SMALL、MEDIUM 和LARGE 三個(gè)數(shù)據(jù)集對(duì)CPU 代碼進(jìn)行了測試，數(shù)據(jù)集的大小根據(jù)網(wǎng)格規(guī)模來區(qū)分（表1）。在三個(gè)數(shù)據(jù)集下，CPU 算法中的主要函數(shù)所占CPU 運(yùn)行時(shí)間比例如表2 所示。在FDTD 的CPU 算法中，大部分時(shí)間消耗在電場向量E 和磁場向量H 修正計(jì)算的多重循環(huán)上。算法執(zhí)行的復(fù)雜度和空間需求估算如下：

　　時(shí)間復(fù)雜度∝ 6×?xí)r間步數(shù)×網(wǎng)格數(shù)×操作數(shù)

　　空間復(fù)雜度∝ 6×網(wǎng)格數(shù)×常數(shù)

　　其中，因數(shù) 6 代表在3D 模型中需要修正的分量的個(gè)數(shù)，常數(shù)表示一個(gè)浮點(diǎn)數(shù)所占的字節(jié)。在串行情況下，由于只能算出一個(gè)點(diǎn)的場值分量，如果要更新整個(gè)網(wǎng)格的場值，就需要一個(gè)三重循環(huán)進(jìn)行計(jì)算。

相關(guān)資料：

基于CPU-GPU 異構(gòu)機(jī)群的FDTD 并行算法加速研究（二）

基于CPU-GPU 異構(gòu)機(jī)群的FDTD 并行算法加速研究（三）