（信号与信息处理专业论文）基于gpu加速fdtd计算速度的研究与仿真.pdf

摘要 时域有限差分方法在电磁计算领域有着广泛的应用。由于f d t d 方法具有天然 的并行特性，它的应用也越来越广泛。然而由于传统的基于c p u 的计算平台在计算 f d t d 方面遇到了瓶颈，c p u 已经无法满足f d t d 计算在速度、内存方面的要求。 图形处理器的出现解决了这一问题。g p u 是专门用来处理并行计算的专用硬件平台， 由于f d t d 的良好的并行特性，正好适合在g p u 上进行计算。 本文针对传统的f d t d 算法的不足，以图形处理器协l ae 1 0 6 0 为计算平台，通 过理论分析和数值模拟，研究并实现了基于g p u 的一维与二维f d t d 算法。同时， 与在c p u 上的一维与二维f d t d 的仿真结果结果进行了比较。得出基于g p u 来进 行f d t d 运算的方法是正确的，并且g p u 确实起到了加速f d t d 计算速度的目的。 同时在划分为不同网格数目下，比较了g p u 与c p u 的计算时间，结果说明仿真数 据越多，g p u 比c p u 的优势越明显。 关键词：时域有限差分法；中央处理器；图形处理器；并行计算 a b s t r a c t t h ef i n i t ed i f f e r e n tt i m ed o m m nm e t h o dh a sb e e nw i d e l yu s e di nt h ef i e l do f t h ec a l c u l a t i o ne l e c t r o m a g n e t i c a st h ef d t dm e t h o dh a st h en a t u r ec h a r a c t e r i s t i co f p a r a l l e l i s m ，s o i t sa p p l i c a t i o nh a sb e c o m em o r ea n dm o r eu n i v e r s a l b u tt h e c a l c u l a t i o nt h a tu s e dt h ec p u si nt h ef d t dh a sm e ts o m ep r o b l e m s t h ec p u sh a v e n ol o n g e rm e e tt h er e q u i r e m e n t so ft h em e m o r ys p a c ea n dt h es p e e do ft h ef d t d c a l c u l a t i o n t h ee m e r g eo ft h eg r a p h i cp r o c e s s i n gu n i t ( g p u ) s o l v e st h ep r o b l e m s g p ui sah a r d w a r et h a ts p e c i a l l yd e s i g n e df o rt h ep a r a l l e lc o m p u t i n g b e c a u s eo f t h eg o o do ft h ef d t di np a r a l l e l i s m , j u s tr i g h tf o rt h eg p u s a i m e da tt h el a c ko ft h et r a d i t i o n a lf d t d a l g o r i t h m ，i nt h i sp a p e r , w ea p p l yt h e g p uc o m p u t i n gp l a t f o r mo ft h et e s l ac 10 6 0 t h r o u g ht h et h e o r e t i c a la n a l y s i sa n d n u m e r i c a ls i m u l a t i o nis t u d i e da n d i m p l e m e n t e d f d t d a l g o r i t h m f o rt h e o n e - d i m e n s i o na n dt w o d i m e n s i o n su s i n gt h et e s l ac 10 6 0 w ec a nc o n c l u d e 廿l a tu s i n g t h eg p ut oc a l c u l a t et h ef d t di sr i g h t i na d d i t i o n ，t h eg p u r e a l l yp l a yt h er o l eo f a c c e l e r a t i n gt h ef d t dc a l c u l a t i o ns p e e d a tt h es a m et i m e ，ic o m p a r et h ec o m p u t i n g t i m eo ft h eg p ua n dt h ec p uu n d e rt h ed i f f e r e n tn u m b e r so ft h eg r i d s t h er e s u l t s s h o wt h a tt h em o r ed a t e sw eh a v et os i m u l a t e ，t h em o r ea d v a n t a g eg p uo v e rt h ec p u k e yw o r d s ：f d t d ；c p u ；g p u ；p a r a l l e lc o m p u t a t i o n 青岛大学硕士毕业论文 1 1 课题背景 第一章引言 在英特尔的奔腾时代，摩尔定律让处理器的主频时钟速度在短短几年显著提 升。随着摩尔定律的进步，从4 5 纳米到3 2 纳米，c p u ( c e n t r a lp r o c e s s i n gu n i t ， 中央处理单元) 集成的晶体管越来越多，能耗越来越大，但在高性能计算方面还 不能满足用户持续增长的需求。同时，基于c p u 的大型计算服务器也由于高额 的采购成本和维护费用成为了无法逃避的弊端。c p u 的设计宗旨就是要从指令 流中得到最高的处理效能，即要用最短的时间完成一项任务，这就是串行计算。 串行计算架构的c p u 天生就不适合高性能计算领域所强调的并行计算工作，这 是为什么c p u 现在成为了高性能计算性能提升的瓶颈。相形之下，始于图形处 理的g p u ( g r a p h i cp r o c e s s i n gu n i t ，图形处理器) 一开始就是被设计来处理大规模 的并行计算，数千万的多边形和像素的并行处理，使得g p u 在发展的同时不断 地加入更多的处理核心并完善并行处理能力。在n v i d i a ( 英伟达) 新的构架下， 可以实现g p u + c p u 的异构计算。两者的协同比以往单纯c p u 运算高出几十倍 甚至几百倍、上千倍。目前，g p u 在市面上已经有几亿之多，而支持 c u d a ( c o m p u t eu n i f i e dd e v i c ea r c h i t e c t u r e ，统一计算设备构架) ，编程的g 8 0 系 列以上的也有l 亿之多。这使大规模并行程序大众化成为可能。自2 0 0 3 年，g p u 多核处理器在浮点性能上开始远远超越c p u 。现在c p u 的性能提升已经放慢， 但是g p u 却仍在急速提升。2 0 0 8 年，多核g p u 与多核c p u 的峰值浮点计算能 力的比率已经超过1 0 。现在最新的g p u 系列t e s l ac 1 0 5 0 的单精度运算峰值为 9 3 3 g f l o p s ，双精度运算峰值为7 8 g f l o p s ，具有2 4 0 个处理器核心；t e s l ac 2 0 5 0 的双精度运算峰值更是达到了5 2 0 6 3 0 g f l o p s ，具有9 6 0 个处理器核心，功耗 只有最新的四核c p u 的1 2 0 ，搭建成本为1 1 0 。迄今为止，如此巨大的计算差 距已经使众多的软件开发商将他们开发的软件中计算密集的部分放在g p u 上来 执行。 f d t d ( f i n i t ed i f f e r e n c et i m ed o m a i n ，时域有限差分) 是在时域内求解麦 克斯韦方程的直接方法。经过几十年的发展，目前已经发展成为一种比较成熟的 数值计算的方法。它的应用范围也越来越广泛，其中有天线的设计、微波电路的 设计、电磁兼容分析、电磁散射计算以及光子学应用等等。无论是从方法本身来 说，还是从具体的研究应用方面来讲，f d t d 目前的发展可以说都是十分的迅速。 在最近的几十年中，关于f d t d 的论文的增长是十分快速的，几乎是以指数的形 式增长。同时，f d t d 与其它的计算有所不同，我们可以把时域有限差分的计算 过程划分为很多个子计算过程，各个子计算过程之间进行着相似的计算。在这些 引言 相似的子计算过程中只存在边界的交换。就这方面而言，f d t d 可以说是一种天 然的并行算法。利用这一点，就可以用g p u 来进行f d t d 计算。从而大大加快 f d t d 的计算速度。不仅如此，g p u 还可以满足f d t d 计算时所需要的极大的 内存，从而能够解决很多在单台p c 条件下根本无法解决的非常复杂的问题。与 巨型机和p c 集群相比，g p u 具有低能耗低成本高性能的优势。 1 2 时域有限差分方法简介 1 2 1时域有限差分的发展 时域有限差分方法是求解m a x w e l l 微分方程的直接时域方法，经过三十多年 的发展已经成为一种成熟的方法，应用范围也越来越广，近几十年来，每年发表 的论文几乎按指数增长。【l 】以下简单的介绍一下f d t d 的发展情况。 l 、y e 9 2 】第一次提了出麦克斯韦方程的差分离散形式，并将其用来进行电磁 脉冲的传播和反射等方面问题的处理。 2 、t a y l o 一3 】用时域有限差分的方法来进行非均匀介质体的电磁散射的分析， 并提出了用吸收边界的方法来吸收外向行波，只是这个吸收边界采用的是比较简 单的插值方法。 3 、m 一4 】提出了一种十分有效的吸收边界条件：在计算区域截断边界处的一 阶和二阶吸收边界条件及其在f d t d 的离散形式，获得了广泛应用。 4 、c h o i 和h o e f e r 5 】用f d t d 分析了波导腔体的谐振问题。 5 、k a s h e r 和g e g 们、m e i 等分别提出了亚网格技术、共形网格技术。 6 、s u l l i v a n 7 用f d t d 研究了生物电磁学的问题，用6 0 - , 7 0 m h z 电磁波照射 下透入到人体内部的电磁场。 7 、b n d 8 】首次给出了时域远场的结果，但论文未给出具体的外推方法； 等提出了三维f d t d 时域近远场外推方法；随后l u e b b 铘提出了二维f d t d 时 域近远场外推方法。【9 】 8 、l a r s o n 1 0 】等进行了能适用于f d t d 的专用计算机方面的研究。从而能够 方便的进行计算电磁波与电大尺寸物体的相互作用。 9 、l u 而b e f s 【1 1 l 研究了用f d t d 的方法处理色散介质的问题。 1 0 、m a l o n e y l l 2 1 进行了在f d t d 中应用阻抗边界条件等方面的研究。 1 1 、s u i t l 3 l 提出用二维f d t d 来计算有集中参数原件的数字和微波模拟电路 模型。 1 2 、b e r e n g e r i l 4 】提出了一种全新的吸收边界。将麦克斯韦方程扩展为场分量 分裂形式，并构成完全匹配层。s a c k s 和g e d n e y 根据b e r e n g e r 的研究又提出了 各向异性介质的完全匹配层，它的支配方程是基于各向异性介质的麦克斯韦方 2 青岛大学硕士毕业论文 程；在f d t d 中这两种吸收边界现在已经得到了广泛的应用。 1 3 、p r a t h e r 和s h i 15 】分析轴对称衍射透镜，给出波长为l 微米平面波和高斯 波入射时，直径为1 0 2 4 7 微米衍射透镜的光波传播特性。 1 6 1 1 2 2f d t d 的应用 f d t d 在电磁研究的多个领域获得广泛的应用【l 】，其中有： 1 、可以用来对辐射天线进行分析：常用的天线有柱状天线，锥形天线，呈 喇叭状的天线，微带天线，使用在手机上的天线，缝隙状的天线，螺旋状的天线 和阵列天线。 2 、微波器件方面以及导行波结构方面的研究：波导和介质波导，微带传播 波导中的孔隙耦合。 3 、散射和雷达截面的计算方面的研究：导体、介质物体以及具有复杂的结 构和形状物体的雷达截面，导弹引导头的电磁波投入分支，我们的身体对电磁波 的吸收，地面下的物体散射等方面。 4 、周期结构分析：频率表面的选择、光栅的传输特性、周期阵列天线和随 机粗糙表面等方面。 5 、进行电子封装、电磁兼容的分析：多线传输和高密度封装时的数字信号 的传输，分析环境结构对元器件性能的影响。 6 、核电磁脉冲的传播和散射，在地面的反射及对电缆传输线的干扰。【1 7 】 7 、光学器件中的光的传播特性和衍射特性。 1 3 论文研究的主要内容 在f d t d 的计算中，计算速度一直是人们所追求的，但是传统的基于c p u 的计算已经遇到了瓶颈，本文介绍了一种基于另一种硬件平台乜p u 的f d t d 计算。 本文首先介绍了f d t d 的发展及应用，f d t d 算法的基本原理等。 其次，介绍了图形处理器的发展历史，对图形处理器在国内外的应用以及 g p u 平台相对与c p u 平台的优势进行了说明，指出了g p u 对于大规模的数据 平行计算的独到之处，但是在处理串行计算时仍需要c p u 的帮助。 最后介绍了g p u 平台的编程模型，并对g p u 加速f d t d 的情况做了仿真， 说明g p u 确实是具有很大的优势。 3 第二章f d t d 的基本理论 第二章f d t d 的基本理论 麦克斯韦( m a x w e l l ) 于1 8 7 9 年建立了电磁场的基本方程，截止到现在为止 电磁波理论以及其应用的发展已经经过了百年的历史。当前，有关电磁波的各个 方面的研究已经在各个领域有了较为深入的研究。此外，它的应用也已经相当的 广泛，比如电波的传播、光纤通信、移动通信、雷达技术、微波、天线、电磁成 像，等等。 电磁波在实际环境中的传播过程十分复杂，例如各个复杂目标的散射，复杂 结构天线的辐射，在波导和微带结构中的传播，实际通信中城市环境、复杂地形 及海面对电磁波传播的影响，等等，具体实际的研究电磁波的特性有着十分重要 的意义。【1 6 】实验及理论分析计算是两个重要的研究手段，二者又是相辅相成的。 其中的分析计算的方法需要用到实际环境中的电磁参数，并用这个参数来求 解麦克斯韦方程边值。但是在通常情况下，解析解仅仅存在于一些经典问题中。 解析解对于计算来说具有重要的指导意义。但是，实际环境是具有复杂性的，往 往需要通过数值解得到具体环境下的电磁波特性。 当今计算机技术的飞速发展，已经为解决麦克斯韦方程提出了许多有意义的 数值方面的方法。例如矩量法( m o m ) 、有限元法( f e m ) 、边界元法( b e m ) 及时域有限差分( f d t d ) 方法【1 8 】等等。并且，随着电磁波应用的广泛和计算机 技术的发展，各种方法的研究也更加深入。 1 9 6 6 年k s 涮1 9 】提出了一种有关电磁场数值计算的新方法一一 f d t d ( f i n i t ed i f f e r e n c et u n ed o m a i n ，时域有限差分) 方法。这种方法是对电磁场 e ，日这两个分量在空间上以及时间上进行交替抽样这种离散方式。每一个e ( 或 者日) 场分量周围有四个日( 或e ) 场分量环绕，应用这种离散方式将含时间变 量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间 电磁场。【2 0 l 这种方法中使用的抽样方式就是我们现在所说的y e o 元胞。时域有限 差分方法可以在时域内直接求解麦克斯韦微分方程。在这个计算过程中，我们把 在空间某个样本点的电场( 或磁场) 与在其周围格点的磁场( 或电场) 进行直接 的关联。与此同时，给每一个元胞都赋予介质参数。我们可以应用这种方法来解 决具有比较复杂的形状的目标的电磁散射、辐射问题。还可以解决非均匀介质物 体的电磁散射、辐射的问题。 2 1f d t d 基本算法 f d t d 方法是一种基于有限差分的方法。这种方法是在时间域和空间域来求 解麦克斯韦方程的两个旋度方程( 2 1 ) 、( 2 2 ) 的电磁分布的数值方法。将安培 方程和法拉第方程这两个方程利用中心极限差分方法进行离散。接下来利用蛙跳 4 青岛大学硕士毕业论文 式的方法在空间领域和时间领域进行交替的计算，通过时间领域的更新来模仿电 磁场的变化，达到数值计算的目的。啪】 v x h ：l _ d + 了 出 ( 2 1 ) v ：一塑+ 了。 a ( 2 2 ) 其中占为电场强度，单位为伏特米；d 为电磁通密度，单位为库伦米；h 为磁场强度，单位为安培米；五为磁通量密度，单位为韦伯米2 ；，为电流密度， 单位为安培米；厶为磁流密度，单位为伏特米2 。各项同性介质中的本构关系 为： 胪占e ，b f f i f h ，j f f i - a e ，吒h ( 2 3 ) 其中s 表示介质介电系数，表示磁导系数，仃表示电导率，表示导磁 率， 仃、o m 分别表示介质的电损耗和磁损耗。真空中盯= o ，吒= 0 ，以及 g = 占o28 8 5 x 1 0 1 2 f m 弘= p o = 4 n 1 0 一h m 直角坐标系下，式( 2 1 ) 、( 2 2 ) 可以写为： 亟一一a l l y ：g 堡+ 盯臣 砂 龙研 “ 孕一孕：占拿+ 呱 ( 2 4 ) 瑟缸所 ， 8 h v o h x a e j 一= 占2 + a 1 5 鲁一鲁= 叫等一q j 一= 一“_ 二一r r 仃 却瑟 。 a f ”5 警一警= 一警一q l j = 一“_ 二一仃仃 貌a x a t m 誓一鲁一，1 1 警一皿1 - 1 o 一o = 一_ 三一仃 a x 执8 t 5l ( 2 5 ) 以下对式( 2 4 ) 、( 2 5 ) 进行f d t d 差分离散。我们假设f ( x ，y ，z ，t ) 表示e 或者日在直角坐标系下的某一分量，在时间域和空间域中的离散取下面的符号 5 第二章f d t d 的基本理论 来表示： f ( x ，y ，z ，t ) = f ( i a x ，j a y ，k a z ，n a t ) = f ”( f ，j ，k ) ( 2 6 ) 对关于时间以及空间的一阶偏导数取中心差分近似，即： 钞( 工，y ，z ，f ) 础 a f ( x , y ，z ，f ) 砂 钞“y ，z ，f ) 砂 舢。竺趋掣1 椭竺型1 掣1 。厂“o ，七+ 三) 一厂”o ，七一三) 一池_ 兰五上 学卜噬掣 ( 2 7 ) 在f d t d 的离散中，电场以及磁场的各节点空间排布如图2 1 所示，这就是 著名的y e e 元胞。由图可见每一个磁场分量由四个电场分量环绕；同样每一个电 场分量由四个磁场分量环绕。【2 l 】这样的电磁场分量的空间环绕与法拉第感应定律 以及安培环路定律的自然结构是相符的。而且非常适合m a x w e l l 方程的差分计 算，能够恰当的描述电磁场的传播特性。此外，电场和磁场在时间顺序上交替抽 样，抽样时间间隔彼此相差半个时间步，应用这种离散方式将含时间变量的麦克 斯韦旋度方程转化为一组差分方程，从而可以在时间上迭代求解，从而不用进行 矩阵求逆运算，因而，由给定的相应电磁问题的初始值，f d t d 方法就可以逐步 地求得以后各个时刻空间电磁场的分布。 2 2 】并在时间轴上逐步推进的求解电磁 场。 在y e e 元胞中的占、日两个分量的空间节点以及时间步取值的整数和半整数 如表2 1 所示。 由上述f d t d 差分方程组可以得到电磁场的时域推进计算方法，如图2 2 所 示。 6 青岛大学硕士毕业论文 图2 1y e e 方案中电场和磁场的位置 表2 1y e e 元胞中e 、h 各分量节点位置 空间分量取样 电磁场分量时间轴t 取样 x 坐标 y 坐标 z 坐标 e l i + l 2 j k e 节点 e p l j + l 2 k n e l j k + l 2 h x l j + l 2 k + l 2 h 邙点 h y i + l 2 j k + l 2n + l ，2 h z i + l 2 j + l 2 k 7 第二章f d t d 的基本理论 f 、h 以及电磁参数初始化 【 l 设置激励源 l 计算空间各处的f 场 j 计算空间各处的，i 场 l 边界处理计算空同各处的f 场 太 一 结果分析 图2 2f d t d 算法的流程示意图 2 2 直角坐标中的f d t d 在这篇论文中我们要讨论的是一维和二维问题，因此这部分要介绍的是如何 构造一维、二维f d t d 在直角坐标下的方程。 2 2 1二维情形 对于二维问题，假设所有物理量均与z 无关，即0 o z 卸，于是式( 2 1 ) 、( 2 2 ) p a y = s 鼍瑚。 一警= 占鲁蚂 ( 波，横电波) 陲一鲁= 叫警一吒皿 8 ( 2 8 ) 青岛大学硕士毕业论文 卢加一盟a t + 以 a z , a , = , a 讲g ) ，+ 鸣 ( t m 波，横磁波)( 2 9 ) b 胁a h , 砂= 噜+ 仃以 由式( 2 8 ) 与( 2 9 ) 我们可以看出，在二维坐标的情形下，电磁场的几个 直角分量可以被划分为两组独立的方程：e 、耳、皿为一组，称为对于e ：的横 电波( t e 波) ；以、以、置为另外的一组，称为对于乞的横磁波( t m 波) 在应用f d t d 方法进行空间离散时，二维坐标f d t d 中曰，h 各分量空间取值 与时间步的取值的整数与半整数如表2 2 所示。电场以及磁场有关节点的空间排 布如图2 3 所示。 表2 2t e 波与t m 波y e e 元胞中层，日各分量节点位置 空间分量取样 电磁场分量时间轴t 取样 x 坐标 y 坐标 h t i + l 2 j + l ：2 n + l 2 t e 波 e 。 i + l 2 j n e y l j + l 2 n 疋 l j n h x l j + l 2 n 十l 2 t m 波 h y i + l 2 j n + l 2 9 第二章f d t d 的基本理论 ( i l j ) + 凹壁! 堑型生盟堕型竺坠些幽 ( 2 1 0 一( m ) 型塑竺韭塑警墨业幽 q jd 啦 2 ) - c q 伽弘l 刚川2 + 长剐+ l 2 i 似化 lyi 矿2 瓴川2 ) = c p ( 册旭叫2 ( f ，川2 ) 一c q ( 小) 盟学( 2 r r l - m ( f + l 2 护c p ( m ) 町( 川2 伊c q ( m ) 盟气芋业( 2 i o 青岛大学硕士毕业论文 坷州2 0 ，_ ，) = 翻( 枘霹o ，) + 四( 神l 盟型攀a x一盟型产 ( 2 ” iyi 式( 2 1 0 ) ( 2 1 5 ) 中： s ( ，1 ) 仃( 册) 1 一墅! a t 黝( 聊) 2 巫a t 砸2 = 砸2 6 ( m ) ( 2 1 6 ) a t 。 2 。2 6 ( m ) c b ( m ) - 硒12 盏 眩 a t2 2 c ( m ) 丛堕一型1 8 = ( m ) 一a t c p ( 朋) 2 巫a t 覃面2 = 触2 9 ( m ) a t2 2 9 ( m ) ( 2 1 8 ) 缸 ：丛碧_ ( 2 1 9 ) l + 盟m ) 篁 2 p ( m ) 以上式( 2 1 0 ) 一( 2 1 5 ) 即是二维坐标情况下的f d t d 电磁场的时间步上 的计算公式。 此外，还可以将t e 波与t m 波在形式上变成更加相似的形式，可将式( 2 1 0 ) 中的相应的标号移动l 2 ，在坐标系上把x 坐标和y 坐标都移动半步，与此同时， t 也移动半步即a t 2 。 利用横电波与横磁波的对应关系，即： r 专，o t t o 一占，巳_ o r ( 2 2 0 ) ej h 。hj e 在编程时就可以将t e 波与t m 波的程序统一起来，二者的程序流程图如图 2 4 所示。 第二章f d t d 的基本理论 ( a ) t m 波 2 2 2 一维情形 ( b l t e 波 图2 4 二维情形f d t i ) 计算电磁场的逐步推进 在一维的f d t d 情况下，可以假设t e m 波是沿z 方向传播的。因此所有的 参数、常量都是与x ，y 这两个参数无关的，即a o 缸= 0 ，a 砂= 0 。在这种情况 下有： 一堡：s 堕+ 仃巨 ( 2 2 1 ) a z a t l 鲁= 叫警一仃q 娩2 2 ， 在一维f d t d 的情形下，电场以及磁场的空间节点如图2 5 所示。一维的 f d t d 的y e e 元胞中的e 、日分量及其时间步的取值的整数和半整数如表2 3 所 示。 1 2 青岛大学硕士毕业论文 ht 一o 图2 5 一维情况场分量空间节点取样 表2 3 一维y 砖元胞中各场分量节点位置 电磁场分量 空间分量取样时间轴取样 e x kn h y 1 c + l 2n + 1 2 式( 2 2 1 ) 和( 2 2 2 ) 的f d t d 离散为： e x n + l ( 栌c a ( 朋) e c a ( 聊) 翌n + l 竺2 坐掣型坠塑 ( 2 2 3 ) z 矿( m 2 ) = 2 c p ( 掰) 矿( 川2 ) 一c q ( m ) 业卷型( 2 2 4 ) 式( 2 2 3 ) 、( 2 2 4 ) 中，c a ( m ) 、c b ( m ) 、c p ( m ) 、c q ( m ) 的含义与式( 2 1 6 ) ( 2 1 9 ) 相同。 2 3数值稳定性 数值稳定性是以有限差分方程来代替麦克斯韦旋度方程，即以差分方程的解 来代替原来的电磁场偏微分方程组的解，只有离散后差分方程的解是收敛和稳定 的，这种代替才有意义。【2 3 】收敛性是当离散间隔趋于零时，差分方程的解在空间 任意一点和任意时刻都一致趋于原方程的解。稳定性是指寻求一种离散间隔所满 足的条件【2 4 1 。在这种条件下，差分方程的解与原方程的解之间的差是有界的。 在本文中我们用二维t e 波来介绍f d t d 方法的稳定性条件。此时的t m 波 的f d t d 公式可写为： 彤川2 ( f ，j + 1 2 ) ：刈2 ( f ，j + 1 2 ) 一祟( 掣( f ，j + 1 ) 一掣( f ，) ) ( 2 2 5 ) “y y n + i 2 ( i + 1 2 , 胪彬圳2 ( i + 1 2 , 卅老( 霹( f + 1 ，舻掣瓴川 ( 2 2 6 ) e + l ( f ，d ：霹“d + 垒壬窆墨兰坚塾笔兰窆燮一旦= i 2 旦卫型垒罢望兰堕型竺 ( 2 2 7 ) 6甑 曲 第二章f d t d 的基本理论 把式( 2 2 5 ) 一( 2 2 7 ) 用时间的本征值表示出来，能够写为： 堡堑韭丝邕幽：ah：(i，+12)at 笪垫坐攀堡塑竺盟：五h y ( i + l 2 , _ ，) a t _ 华= 掣2 瓴力 用v 代表电场与磁场分量，然后能够得到： 竺：! ! ：二匕：! ! ：五y 。 a t ( 2 2 8 ) ( 2 2 9 ) ( 2 3 0 ) ( 2 3 1 ) q = v - + 2 为增长因子，带入式( 2 2 4 ) ，然后方程两边同时除以y ”抛， 可以得到方程： q 2 一勉幻一l = 0 ( 2 3 2 ) 式( 2 3 2 ) 的解为： g = t 2 a t 小+ ( 等) ( 2 3 3 ) 因为要求电场不发散，所以i q i 1 ，要满足该条件，只要： r c ( 五) = ”2 f si m ( 五) 云( 2 3 4 ) 时间本征值必须在数轴的负半轴的范围内才能满足其稳定性的条件。 考虑在均匀介质中传播的电磁波，方程解得平面波形式为： i e ( f ，) = 毛p 肌恤+ b 抄础 ( 2 3 5 ) l 彤( f ，_ ，) = h 知少心黼抄删 ( 2 3 6 ) 【彬( f ，_ ，) = 民e j o 黼抄础 ( 2 3 7 ) 将式( 2 3 5 ) - - ( 2 3 7 ) 带入式( 2 2 5 ) 一( 2 3 1 ) 可以得到f d t d 方法中的特征 值： a 2 = 一：兰 警+ 警 ( 2 3 8 ， 式( 2 3 8 ) 中兄为纯虚数，可以看出对所有t 、七，、a 满足下列条件： 1 4 青岛大学硕士毕业论文 r e c a ，= 。，之v 厮s - m c a ，2 v 膈 c 2 3 9 ， 式中v = l i ，为光在介质中传播的速度。 利用式( 2 3 8 ) 和( 2 3 9 ) 可以得到： a ts 式( 2 4 0 ) 称之为二维情况下f d t d 的c o u r a n t 稳定条件。 2 4f d t d 中常见的激励源 ( 2 4 0 ) 在应用f d t d 的方法进行电磁问题分析时，要面对激励源的选择的问题。就 是说需要选择一个合适的入射波，然后将入射波变为适当的形式迭代到f d t d 中。激励源可以分为随时间周期性变化的时谐场源，脉冲形式的场源；面源、线 源、点源等。下面介绍本文仿真所用到得激励源：高斯脉冲源和时谐场源。 2 4 1时谐场源 黔n 。刎 这个场源是一个自t = 0 开始的半无限正弦波列。 2 4 2 高斯脉冲源 高斯脉冲的时域形式为： 眦) = e x p ( 等竽) ( 2 4 1 ) ( 2 4 2 ) 高斯脉冲的峰值出现的时刻为t = t o 时刻，它宽度由f 这个参数来决定。上面 的式子的傅里叶变化为： w ) = 三e x p ( - 伽死一竿) ( 2 4 3 ) 2 5 激励源的设置 一般情况下主要用下面这两种方式来为f d t d 引入激励源： 1 5 第二章f d t d 的基本理论 1 、激励源替代法： 激励源替代法是采用源处的场分量去替代y e e 网格计算的场分量。对于t e 波而言，它的电场只有巨这个分量。在( ，工) 点采用频率为的正弦脉冲，e 按以下形式变化： f ( ，五) = e os i n ( 2 7 r f o n a t ) ( 2 4 4 ) 激励源替代法由于形式简单，可以很方便的使用，但是缺点是很容易产生虚 假反射，不利于计算。 2 、激励源叠加法 激励源叠加的方法与替代法的不同之处在于不是替代激励源而是将激励源 叠加在y e e 网格所计算的场分量上。 例如对于t e 波而言，在( ，五) 这一点采用频率为厶的连续的正弦波，e 的 f d t d 差分方程为： 哎。( ，工) = e os i n ( 2 a f o n a t ) ( 2 4 5 ) 胪眺+ 譬( 型坠堕罢盥型趔 一型塾址塑警迅幽) + 刚n 似- ) ( 2 4 6 ) v 采用叠加激励源的方式，可以消除虚假反射的问题。 2 6 吸收边界条件 由于现在的计算机容量和计算机计算能力都是有限的的，f d t d 的计算目前 只能在有限区域进行。因此，在计算的时候必须把f d t d 的计算空间截断。为了 能够模拟在开放区域内的电磁散射与辐射过程，在f d t d 计算区域的截断边界处 吸收边界条件必须要清楚，并且这种截断边界在截断处不能产生比较明显的电磁 波反射。吸收边界条件一开始采用的是简单的插值边界。直到1 9 8 1 年，g m u r 提出了m u r 近似吸收边界条件【2 4 】之后，f d t d 才能够开始解决实际的问题。1 9 9 4 年b e r e n g e r 提出了二维理想匹配层( p e r f e c t l ym a t c h e dl a y e r ，p m l ) 吸收边界 条件【2 6 】。其能够吸收不同方向不同频率的电磁波。1 9 9 6 年，g e n d e y 提出了各向 异性介质p m l 理论【2 刀。与m u r 二阶吸收边界条件比，p m l 可以说是目前最好 的吸收边界条件。 2 6 1b e r e n g e r 完全匹配层 1 6 青岛大学硕士毕业论文 完全匹配层( p e r f e c tm a t c h e dl a y e r ，p m l ) 的概念是由b e r e g e r 在1 9 9 4 年 提出得。这种吸收方法的原理是阻止入射波穿过边界进入p m l 层。这就需要在 f d t d 区域截断边界处设置一种特殊的介质层。说它特殊是因为它的波阻抗和与 它相邻的波阻抗是完全匹配的。并且有能够使进入p m l 层的透射波迅速衰减这 种特性。具有这种特性之后，有限厚度的p m l 层也能够具有很不错的吸收入射 波的作用。 以二维t e 波为例，t e 波中只有皿、巨、e ，分量，在直角坐标系中，麦克 斯韦方程写为： 8 e 8 h i 言2 畜 岛鲁一警 眩4 7 ， 陲一鲁一鳓警 在p m l 中，b e r e n g e r 采用的方法是假设将磁场分量h ：分为两个子分量：h 珏 和h 巧，且h z _ h h 巧。这样就可以将将麦克斯韦方程改写为以下形式： 气鲁+ q 乓= 半 岛鲁蚋2 一鼍掣 他4 8 ， 风百0 h = + 也= 一誓 百a l l y + 考 其中吒、q 、分别代表介质的电导率和磁导率。上式中当 吒= = 仃，= = o 时，式( 2 4 7 ) 变为自由空间的麦克斯韦方程式( 2 4 8 ) 。 在这里只对t e 波进行研究，对于t m 波的情形同样适用。 由于在这层p m l 介质中，电磁波会以很快的速度进行衰减，所以在一般情 况下就采用指数差分的方法对上面的式子来进行离散。但是在我们的实际应用 中，指数差分法与中心差分法就二者的精度方面来说是差不太多的，此外二者还 有大体上相当的计算复杂度。 2 6 2 完全匹配层的结构 完全匹配层的基本结构如图2 6 所示。空白部分为计算区域，灰色部分为 1 7 第二章f d t d 的基本理论 p m l 层。在计算区域中，应用常规的计算法计算麦克斯韦方程即可。在计算区 域中的产生的外形波在穿过p m l 层的时候被p m l 层所吸收。 o x p m l ( 0 ，0 ，o i y o 。) 图2 6 二维p m l 边界参数设置 另外p m l 层不能是无限厚度的，只能是有限厚度，通常设置p m l 层数为 l o 层。 2 6 3 二维m u r 吸收边界条件的f d t d 形式 对于二维电磁场r o 题，m u r 之处二阶近似吸收边界条件可降低为只含有五， 分量的一节导数，从而可以令数值计算简单化。对于t m 波而言有： c 三盖专嘉畦弘l - 0 旺4 9 ， 将式( 2 8 ) 的第二个方程带入式( 2 4 9 ) ，变为 l 堕一! 婆一c # 0 2 n xl ：o ( 2 5 0 ) lo x o t co t 22o y o tj ，o 将式( 2 5 0 ) 式对t 进行积分，并设初始时刻场为零，可得： t 警一三鲁一等等k = o ( 2 州， 对于t e 波，有： c 三去专嘉畦等，4 卜 眨5 2 ， 将式( 2 9 ) 式代入，e 式变为： 1 8 青岛大学硕士毕业论文 l 姿一! 粤+ 竺姿l ：o ( 2 5 3 ) 【- 缸a c 西2 2 砂新i x - 0 将式( 2 5 3 ) 式对t 进行积分，并设场的初始值为零，可得 【警一1 cl 西_ 王h 一詈等k 。= o ( 2 5 4 ) m u r 二阶吸收边界条件与一阶近似吸收边界条件相比，m u r 仅比一阶近似多 了一项：一阶导数项，即式( 2 5 3 ) 与( 2 5 4 ) 两式括号中的第三项。它的形式 要简单的多。此外，t e 波和t m 波的二阶吸收边界条件式( 2 5 1 ) 与( 2 5 4 ) 有 对偶关系。 下面讨论m u r 吸收边界条件的f d t d 离散式。先看t m 波时的一阶近似。 由( 吴一! 昙) 厶i ：o ，可得 i 堕一! 盟l ：o l 叙 caj ，1 0 将式( 2 5 5 ) 在( i + 1 2 ，j ) 点和咐1 2 ) a t 时刻做离散，各项为： 曼生r 2 ：垡：塑! ! 盈二垡：塑! 盟 融l i + 1 1 2 , j 6 望r 2 ：垡：堑! ! 兰! 盈二竺堡! ! 兰! 旦 a t l + v z , j a t ( 2 5 5 ) ( 2 5 6 ) 捎去式( 2 5 6 ) 中e 在( i + 1 2 ) 和( n + i 2 ) 半整数点处的值，最后可得到 e + 1 ( f “，) + 竺等【霹+ 1 0 “，_ ，) 一f ( f ，_ ，) 】 ( 2 5 7 ) 对于t m 波m u r 吸收边界条件的二阶形式，其f d t d 的离散式比式( 2 5 7 ) 多一项，即罢二r 胆：旦= 玉塑旦竺生巡：! ! ：坠! 丝塑 ( 2 5 8 ) o y j ，+ l ，2 ， a t t e 波的吸收边界条件的一阶近似于t m 波形式完全相同，即 彤+ 1 ( “) = 彤( “l ，_ ，) + 竺等 彤+ 1 ( l ，_ ，) 一彤( “) 】 ( 2 5 ” t e 波吸收边界条件的二阶近似为： h 二n + l ( = 彤( f + + 瓦c a t l - 8 仃：+ 1 ( 舻川( 例+ 丽c 2 e a t 【e + 坨( f ，j + l 2 ) - e + 1 几( f ，j - l 2 ) + e 州2 ( i + 1 2 ，j + l 2 ) ( 2 6 0 ) 一e + 1 胆( f + l ，_ ，一l 2 ) 1 9 第二章f d t d 的基本理论 对于二维f d t d 的情况，电磁场除了具有e 外，还具有以、日，这两个分 量。在应用时域有限差分方法来计算边界处的t m 元胞彤+ 1 儿( f ，j + 1 2 ) 以及 日儿( i + 1 2 ，力时，并没有涉及到截断边界以外的e 节点或者是日的节点。只 有e + 1 ( f ，) 涉及到了截断边界外侧的日节点。由此可知，仅仅要给出边界处的 切向场分量巨的吸收边界条件就可以了。类似的情况，对于t e 波而言，仅仅需 要给出边界处切向场分量日的吸收边界条件就可以了。 最后，如果所用的计算区域为矩形0 z a ，0 x b ，就会存在有四个截断 边界。以上情况下，我们只讨论了x = 0 的左侧界面的吸收边界条件。另外的三边 具有与上述情形类似的结果。 2 7本章小结 本章系统介绍了f d t d 方法的原理、基本算法的相关问题。这些算法有f d t d 计算的稳定性，吸收边界条件，激励源等，这是进行f d t d 应用的基础及关键。 青岛大学硕士毕业论文 第三章图形处理器g p u 图形处理器( g r a p h i cp r o c e s s i n gu n i t ，g p u ) 是一块专用集成电路芯片是个 人计算机中显示卡的核心部件。负责计算机中有关图像图形显示的场景计算、变 换和输出。n v i d i a 公司在1 9 9 9 年发布g e f o r c e 2 5 6 图形处理芯片时首先提出 g p u 的概念。g p u 使显卡减少了对c p u 的依赖，并进行部分原本c p u 的工作， 尤其是在3 d 图形处理时。 3 1g p u 概述 在二十世纪八十年代的初期，出现了以g e ( g e o m e t r ye n g i n e ) 为标志的图形 处理器。从此，计算机图形学的进展进入了由图形处理器引导其发展的年代。 g e 处理器的核心部分是四位向量的浮点运算。它的一个寄存器定制码可以定制 出各种不同的功能，分别应用于在图形渲染流水线中，可以实现矩阵、裁剪、投 影等各种运算。如果有十二个这种的单元，我们还能够很完整地实现三维图形的 流水线的功能。图形界面g l 及其后的o p e n g l ，成为图形用户界面工业标准。 在二十世纪八十年代和九十年代之前，图形处理还是主要的在c p u 上来完 成。但是，随着g p u 的功能的不断增强以及结构的不断的完善，图形处理便转 移到了g p u 上。我们可以把图形处理的流水线功能划分为下面两个方面的内容。 第一个方面是对图形进行几何变化，还有对图元的属性进行处理；第二个方面主 要是扫描转换进行光栅化以后完成一系列的图形绘制处理。这方面主要包含有各 种光照效果和合成、纹理映射、遮挡处理、反混淆处理。【2 8 】 进入二十世纪的年代，英伟达公司打入个人电脑三维市场，它在一九