CN102207987A - 基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法 - Google Patents

Abstract

本发明提出一种基于OpenCL的GPU加速三维时域有限差分(FDTD)电磁场仿真的方法，应用于电磁场仿真分析领域，利用OpenCL简化FDTD方法对电磁场仿真的设计与实现。首先设置FDTD仿真参数并动态申请内存空间，构建三维的电磁模型并进行网格剖分，然后初始化OpenCL的计算参数，对三维电磁模型基于OpenCL进行FDTD加速仿真，最后用OpenCL标准库函数释放显存中存储的数据，释放CPU内存资源。本发明方法显著提升了FDTD电磁场仿真速度，与利用CPU计算相比速度提升可达5-10倍，且具有UPML吸收边界条件，可以模拟电磁波在自由空间的传播，使得FDTD电磁场仿真具有更广泛的实际应用。

Description

基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法

技术领域

本发明涉及电磁场仿真分析领域，具体涉及一种基于OpenCL(Open Computing Language，开放计算语言)的GPU(Graphic Processing Unit，图形处理器)加速三维时域有限差分电磁场仿真的方法。

背景技术

FDTD(Finite-Difference Time-Domain，时域有限差分)是电磁场仿真计算领域的一种常用方法。FDTD方法是把麦克斯韦(Maxwell)方程组在时间和空间上进行差分离散，通过蛙跳式(Leap frog algorithm)迭代，使空间内的电场和磁场进行交替计算；通过时间上更新来模仿电磁场随时间的变化，实现数值计算目的。用FDTD方法分析电磁学问题时，要事先构建研究对象的空间几何参数、材料参数模型，并可对计算精度，计算复杂度，计算稳定性进行控制。其优点是能够直接模拟空间电磁场随时间的变化和分布情况，计算精度高，是目前应用较多的电磁学数值模拟方法之一。

FDTD电磁场仿真分析方法已被成功应用在电磁兼容分析、天线设计、射频器件设计、微纳光学等研究领域。然而由于FDTD为时域全波电磁学算法，因此面临仿真时间过长的关键性制约。目前解决FDTD仿真时间长的方法主要有利用高性能并行计算机进行仿真计算，然而高性能并行机的成本高昂，因此该方法并没有得到普及。由于GPU具有CPU所不能比拟的高性能并行计算能力，并且成本低廉，因此利用GPU加速FDTD仿真速度成为一种低成本高性能的仿真方法。

利用GPU加速FDTD仿真虽然具有高性能低成本的优点，但需要设计开发人员具备深厚的图形编程理论知识储备。近些年Nvidia公司和AMD公司相继开发了CUDA和Brook+等通用计算程序接口用以推广其GPU应用于科学计算，使设计开发人员不必精通图形编程理论，也可以利用GPU进行科学计算。目前国外已经有学者成功实现了利用上述接口实现GPU加速二维或三维FDTD仿真的例子。但上述FDTD加速方法只能运行在GPU上，不能运行在CPU上，并且上述的FDTD加速仿真实例不具有吸收边界条件，而不具备吸收边界条件的FDTD方法实用性很有限，只能仿真金属波导结构。

发明内容

本发明针对目前FDTD电磁场仿真计算时间过长的弊端，以及目前FDTD仿真加速方法只能运行在GPU上且实用性有限的缺陷，提出了一种基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法。

本发明是一种基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，具体包括以下步骤：

步骤1：设置FDTD仿真参数，为所述仿真参数动态申请内存空间；

步骤2：构建三维的电磁模型，对该电磁模型进行网格剖分；

步骤3：初始化OpenCL的计算参数，包括：获得FDTD仿真空间的大小、设置仿真所用处理器类型、创建上下文、获得工作组的维数及大小信息、创建指令序列及数据缓存区、装载FDTD加速程序文件、声明磁场仿真计算的内核函数以及电场仿真计算内核函数；

步骤4：对三维电磁模型基于OpenCL进行FDTD加速仿真；

步骤5：利用OpenCL标准库函数clReleaseMemObject()释放显存中存储的数据；

步骤6：仿真完成后，释放CPU内存资源。

所述的步骤3中，通过设置仿真所用处理器类型为GPU或者CPU来设置本发明方法是基于GPU的FDTD电磁场仿真或者基于CPU的FDTD电磁场仿真。

所述的步骤4具体通过下面过程进行FDTD加速仿真：设置时间迭代步数；为每一个迭代步设置其对应的开始时间；按照步骤1中设置的FDTD仿真参数设置激励源；执行基于OpenCL的磁场仿真计算内核函数；执行基于OpenCL的电场仿真计算内核函数；保存需要的电磁场数据。

本发明的优点和积极效果在于：(1)本发明的加速方法显著提升了FDTD电磁场仿真速度，与利用CPU计算相比速度提升可达5-10倍；(2)本发明的加速方法具有UPML吸收边界条件，可以模拟电磁波在自由空间的传播，使得FDTD电磁场仿真具有更广泛的实际应用；(3)本发明的加速方法既可以在GPU运行，也可以在CPU运行，并可利用多核CPU的并行性，采用本发明的方法进行电磁场仿真计算，可使其具有跨硬件平台的通用性和可移植性，使得FDTD电磁场仿真得到更大程度的普及。

附图说明

图1是本发明加速三维FDTD电磁场仿真方法的整体步骤流程图；

图2是本发明方法中UPML设置的部分示意图；

图3是本发明的GPU工作组示意图；

图4是基于OpenCL的磁场计算时y轴方向E_x和E_z数据设置的流程图；

图5是基于OpenCL的磁场计算时z轴方向E_x和E_y数据设置的流程图；

图6是基于OpenCL的磁场计算时x轴方向E_y和E_z数据设置的流程图；

图7是基于OpenCL的电场计算时y轴方向H_x和H_z数据设置的流程图；

图8是基于OpenCL的电场计算时z轴方向H_x和H_y数据设置的流程图；

图9是基于OpenCL的电场计算时x轴方向H_y和H_z数据设置的流程图；

图10是平面形式存储电磁场示意图；

图11是本发明实施例中固定迭代时间步数1000，改变FDTD计算空间条件下GPU和CPU计算时间对比图，(a)是GPU和CPU计算时间对比图，(b)是GPU和CPU计算时间比值的曲线示意图；

图12是本发明实施例中固定FDTD仿真空间为120×120×120网格，改变迭代时间步数条件下GPU和CPU计算时间对比图，(a)是GPU和CPU计算时间对比图，(b)是GPU和CPU计算时间比值的曲线示意图；

图13是本发明实施例中固定迭代时间步数1000，改变FDTD仿真空间条件下基于OpenCL的双核并行CPU计算时间与普通CPU计算时间的对比图；(a)是基于OpenCL的双核并行CPU计算时间与普通CPU计算时间的对比图，(b)是普通CPU计算时间与基于OpenCL的双核并行CPU计算时间比值的曲线示意图。

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

OpenCL是一种面向异构平台系统通用目的的并行编程开放式语言标准，利用OpenCL可以简化FDTD方法对电磁场仿真的设计与实现，使GPU加速FDTD方法进行仿真这一技术更加通用。本发明的一种基于OpenCL的GPU加速三维FDTD电磁场仿真的方法，利用GPU并行处理特性并结合OpenCL接口功能实现了具有单轴各项异性完全匹配层(UPML)吸收边界条件的三维FDTD方法的高性能加速仿真。采用本发明的方法进行电磁场仿真具有跨硬件平台的通用性和可移植性，基于OpenCL编译的语言程序可以运行在CPU或GPU硬件上，并可充分发挥多核CPU的并行计算能力；另外本发明方法采用了UPML吸收边界条件，可以很好的吸收传播到边界面处的电磁波，使电磁波不会发生虚拟反射，从而模拟电磁波在自由空间的传播，在根本上扩展了三维FDTD方法的应用范围。

如图1所示，本发明基于OpenCL的GPU加速三维FDTD电磁场仿真的方法，包括下列步骤：

步骤1：参数设置及内存空间动态申请。在进行FDTD电磁场仿真之前需要对电磁场仿真的参数进行设置及其内存空间动态申请。

步骤1-1：FDTD电磁场仿真参数设置，包括设置空间步长、时间步长、UPML吸收边界以及激励源。

步骤1-1-1：设置空间步长Δα。

为了减小FDTD的数值色散，FDTD仿真的空间步长的设置应满足的条件如(1)式。

$\mathrm{\Δ\α}=\frac{{\mathrm{\λ}}_m}{10},$ α＝x，y，z    (1)

式(1)中，λ_m为激励源的最高频率所对应的波长。根据三维空间步长Δα和三维FDTD空间上三个坐标方向的总尺寸L_x，L_y和L_z，确定FDTD空间中在三个坐标方向上的总网格数N_x，N_y和N_z，如式(2)，当不能整除时采取四舍五入方法。

$N_{\mathrm{\α}}\≅\frac{L_{\mathrm{\α}}}{\mathrm{\Δ\α}},$ α＝x，y，z    (2)

步骤1-1-2：设置时间步长。

为了满足FDTD仿真的数值稳定性，时间步长Δt的设置需要满足的条件如(3)式。

$\mathrm{\Δt}=\frac{0.9}{c\sqrt{\frac{1}{\mathrm{\Δ}{x}^2}+\frac{1}{\mathrm{\Δ}{y}^2}+\frac{1}{\mathrm{\Δ}{z}^2}}}---\left(3\right)$

式(3)中c为真空中的光速。Δx、Δy、Δz分别是在三维方向的空间步长。

步骤1-1-3：设置UPML吸收边界。

设置UPML吸收边界在x-y-z三维方向的网格层数n_x-PML1，n_x-PML2，n_y-PML1，n_y-PML2，n_z-PML1和n_z-PML2。本发明实例中，取三维方向的UPML网格层数为10层。一般情况下，UPML吸收边界区域中的参数σ和κ的取值是逐层渐变的。以x方向的参数设置为例，x方向的参数σ_x和κ_x由x＝0到x＝d按照式(4)-(6)计算：

${\mathrm{\σ}}_x\left(x\right)={\left(\frac{x}{d}\right)}^m{\mathrm{\σ}}_{x,\max }---\left(4\right)$

κ_x(x)＝(g^1/Δx)^x    (5)

${\mathrm{\σ}}_{x,\max }=-\frac{(m+1)\ln \left[R\left(0\right)\right]}{2\mathrm{\ηd}}---\left(6\right)$

其中，d为UPML吸收边界厚度，d＝Δxn_x-PML1，Δx为x方向的空间步长，η为波阻抗，R(0)为可事先设定的外行波垂直入射到分界面时的反射系数，参数m通常可取3≤m≤4，参数g通常可取2≤g≤3。其它方向的σ和κ的取值与此类似。除了一般的情形外，还必须单独考虑下面列出的FDTD仿真区域中σ和κ的特殊取值情况：

(1)无耗、各向同性的内部区域：σ_x＝σ_y＝σ_z＝0，κ_x＝κ_y＝κ_z＝1。

(2)UPML吸收边界的6个面区域：

垂直于x方向的两个截断边界面：σ_y＝σ_z＝0，κ_y＝κ_z＝1；

垂直于y方向的两个截断边界面：σ_x＝σ_z＝0，κ_x＝κ_z＝1；

垂直于z方向的两个截断边界面：σ_x＝σ_y＝0，κ_x＝κ_y＝1。

(3)UPML吸收边界的12条棱区域：

平行于x轴的四条棱：σ_x＝0，κ_x＝1；

平行于y轴的四条棱：σ_y＝0，κ_y＝1；

平行于z轴的四条棱：σ_z＝0，κ_z＝1。

(4)UPML吸收边界的8个角区域：σ和κ按照一般情况取值。

步骤1-1-4：激励源设置方案，常用的FDTD仿真空间激励源可设置为电流源与偶极子源。具体的设置方案如下。

电流源设置方法：在麦克斯韦磁场旋度方程的右侧添加电流密度项，如式(7)所示。其中J表示电流密度，

为电流密度的矢量形式。在FDTD仿真中，可根据实际情况，选择在x，y，z三维的某一方向上施加该激励源。

$\▿\×\stackrel{\→}{H}=\frac{\∂\stackrel{\→}{D}}{\∂t}+\stackrel{\→}{J}---\left(7\right)$

其中，

表示旋度符号，

表示磁场矢量，表示对磁场矢量取旋度，

表示电位移矢量，

表示电流密度矢量。电流密度J的时域波形可以是下式(8)-(10)中的一种。

余弦波源：

高斯波源：

$J(i_s,j_s,k_s)=J_0\exp (-{\left(\frac{t-t_0}{t_w}\right)}^2)---\left(9\right)$

高斯余弦波源：

其中，ω₀表示角频率，

表示初始时刻的角度位置，t_w表示脉冲的宽度，t₀表示脉冲中心所在的时间，i_s、j_s、k_s分别表示电流源所在空间位置的坐标，J₀表示电流密度峰值。

偶极子源设置方法：在麦克斯韦磁场旋度方程的右侧施加电偶极子项，如式(11)。在FDTD计算中，电偶极子通常限制在一个立方体网格的体积内，即Δα³。在FDTD计算中，可根据实际情况，选则在x，y，z三维的某一方向上施加该激励源。

$\▿\×\stackrel{\→}{H}=\mathrm{\ϵ}\frac{\∂\stackrel{\→}{E}}{\∂t}+\frac{1}{{\mathrm{\Δ\α}}^3}\frac{d\stackrel{\→}{p}}{\mathrm{dt}}---\left(11\right)$

其中，

表示偶极子源矢量，

表示电场矢量，电场E的迭代公式如(12)所示：

$E^{n+1}=E^n+\frac{\mathrm{\Δt}}{\mathrm{\ϵ}}{[\▿\×\stackrel{\→}{H}]}^{n+1/2}-\frac{\mathrm{\Δt}}{{\mathrm{\ϵ\Δ\α}}^3}{\left[\frac{\mathrm{dp}}{\mathrm{dt}}\right]}^{n+1/2}---\left(12\right)$

ε表示介电常数，n表示第n时间步，电偶极子源p可以是下式(13)-(15)中的一种。

余弦波源：

高斯波源：

$p(i_s,j_s,k_s)=p_0\exp (-{\left(\frac{t-t_0}{t_w}\right)}^2)---\left(14\right)$

高斯余弦波源：

其中p₀表示偶极子源的峰值，i_s、j_s、k_s分别表示电偶极子源所在空间位置的坐标。

步骤1-2：内存空间动态申请。在CPU内存中声明FDTD仿真计算所用到的所有常量和变量，并为其开辟相应存储空间。

内存空间采用动态分配方式，并且数据连续存储。对于三维电磁场数据变量如表示FDTD仿真空间的在x方向、y方向、z方向的电场E_x、E_y、E_z，表示FDTD仿真空间的在x方向、y方向、z方向的磁场H_x、H_y、H_z，以及由UPML引入的三维临时电磁参量，如步骤4中引入的B_x，B_y，B_z，D_x，D_y，D_z，F_x，F_y，F_z等，访问寻址方式为：

φ[i×N_yz+j×N_z+k]，N_yz＝N_y×N_z    (16)

其中i，j和k为分别代表三维坐标系x轴、y轴和z轴的所对应的空间位置点。N_x，N_y和N_z为x轴、y轴和z轴的所对应的网格剖分数。利用此种数据存储方式可以提高数据寻址和读取的效率。

步骤2：三维模型构建及网格剖分：FDTD仿真分析模型的建立可以采用直接在程序中建模或软件自动建模。三维模型的构建及其网格剖分具体包括下面两个步骤。

步骤2-1：根据所要建立的电磁模型的复杂程度，判断模型建立的方式。对于简单电磁模型，可直接进行三维几何模型的建立与网格剖分，然后执行步骤3。

步骤2-2：复杂电磁模型的构建及网格自动剖分。对于复杂的电磁模型，需要借助AutoCAD等建模软件进行三维几何模型的构建，然后再对构建的模型进行网格剖分和电磁材料参数包括介电常数ε、电导率σ和磁导率μ的赋值。

利用Auto CAD构建需要仿真的三维几何模型，或读入已经用CAD三维建模软件建立的三维几何模型，不同电磁材料部件用图层来区分，然后利用Object ARX二次开发工具开发自动剖分程序输出有效网格位置及其电磁材料参数，生成“几何-电磁参数”描述文件，然后FDTD仿真方法调用该文件，将剖分网格和电磁参数读入计算内存，进行FDTD计算。

所述的自动剖分程序的建立方法是：顺次用FDTD仿真空间中的每一个边长与空间步长相等的立方体单元去与所建立的三维几何模型求交集。如果该立方体单元在目标内部则相交体为立方体单元本身；如果在边缘则相交体为立方体单元的一部分，如果在外部则相交体为自由空间。据此，当相交体的体积不小于立方体的一半时，则可以认为该立方体单元被目标材料完全填充，反之则认为该网格为自由空间。在“几何-电磁参数”描述文件中，可以用序号“0”来表示自由空间，用序号“1”、“2”、“3”、…等来表示不同部件材料。

步骤3：初始化OpenCL计算参数。

步骤3-1：获得基于GPU的FDTD仿真空间的大小。与CPU内存的存储方式不同，GPU的存储单元允许float4型数据结构的声明，即允许创建一种包含4个float型数据的新的数据类型。因此利用float4型数据，可以将E_x，E_y，E_z，H_x，H_y，H_z等数组沿Z轴以每4个相邻数据为一组，从而将数据量压缩4倍。因此基于GPU仿真的总的FDTD仿真空间的大小为N_x×N_y×N_z/4。

步骤3-2：设置仿真所用处理器类型。因为基于OpenCL标准库的程序既可以运行在CPU上，又可在GPU上运行，因此，需要在OpenCL初始化阶段，判断所用处理器类型。OpenCL标准库中CL_DEVICE_TYPE_CPU和CL_DEVICE_TYPE_GPU分别代表仿真计算所用的设备为CPU和GPU，本发明方法主要针对利用GPU加速计算FDTD，在这里将仿真所用处理器类型设置为CL_DEVICE_TYPE_GPU。若是在CPU上运行，则设置仿真所用处理器类型为CL_DEVICE_TYPE_CPU，并可利用多核CPU进行并行计算。

步骤3-3：创建OpenCL上下文。利用OpenCL标准库函数中的函数clCreateContextFromType()创建上下文，并利用函数clGetContextInfo()获得上下文信息。

步骤3-4：检查OpenCL设备。利用OpenCL标准库函数clGetDeviceInfo()获得OpenCL的工作组(workgroup)的维数及大小等信息。

步骤3-5：创建OpenCL指令序列。利用OpenCL标准库函数中的clCreateCommandQueue()创建指令序列，并创建E_x，E_y，E_z，H_x，H_y，H_z等三维电磁场数据变量以及由UPML引入的三维临时电磁参量的数据缓冲区。

步骤3-6：创建OpenCL程序对象。OpenCL的FDTD加速程序文件以“cl”为后缀，如“＊.cl”。运行“＊.cl”文件将执行FDTD加速计算。利用OpenCL标准库函数clCreateProgramWithSource()将基于OpenCL的FDTD加速程序文件装载，并利用函数clBuildProgram()将加速程序文件转化为GPU内核可执行的程序文件。

步骤3-7：创建声明基于OpenCL的磁场仿真计算的内核函数。所述的磁场仿真计算的内核函数保存在“＊.cl”文件中，其具体实现了步骤4-4中的内容。

步骤3-8：创建声明基于OpenCL的电场仿真计算的内核函数。所述的电场仿真计算的内核函数保存在“＊.cl”文件中，其具体实现了步骤4-5中的内容。

步骤4：进行基于OpenCL的FDTD加速仿真。

步骤4-1：设置时间迭代步数n，则总的迭代时间为t＝n×Δt。

步骤4-2：开始循环迭代，设置每一迭代步所对应的开始时间为t₀+Δt×(n-1)。其中t₀为计算起始时刻，一般情况t₀＝0。

步骤4-3：设置激励源，按照步骤1-1-4中的设置方案进行设置。

步骤4-4：执行步骤3-7中声明的基于OpenCL的磁场仿真计算内核函数，具体为：

步骤4-4-1：在时间步n＝1时，设置FDTD仿真所用的内核参数，利用OpenCL标准库函数中的clSetKernelArg()函数取得FDTD磁场计算所用三维电磁场数据变量以及三维临时电磁参量的存储地址。

步骤4-4-2：当利用GPU进行计算时，利用OpenCL标准库函数clEnqueueWriteBuffer()，将指令序列数据缓冲区中FDTD磁场计算所用的三维电磁场数据变量以及三维临时电磁参量等传递到GPU显存。

步骤4-4-3：利用OpenCL标准库函数clEnqueueNDRangeKernel()激活GPU内核进行FDTD磁场加速仿真。基于OpenCL内核函数的FDTD磁场加速仿真如下步骤。

步骤4-4-3-1：OpenCL标准对于GPU计算引入了工作组(Work Group)和本地变量(Local)的概念，其中工作组是由若干工作项目(work item)组成的一片计算空间。工作组和本地变量是相对应的概念，工作组中的工作项目(Work Item)直接利用本地存储空间(Local Memory)进行数据的存储操作和计算操作。在同一工作组中的工作项目共享一片本地存储空间，并可执行并行计算。

对于OpenCL标准来说，在存储的类型上分为本地存储空间(1ocal memory)和全局存储空间(global memory)。全局存储空间是对于所有工作组和工作项目可见的。相对于全局存储空间来说，本地存储空间是指在同一工作组内可见的。同一工作组内的工作项目可以共用一片本地存储空间，本地存储空间可映射到全局存储空间。

本发明方法中工作组大小设置为n_x×n_y×n_z＝2×2×30。

步骤4-4-3-2：定义全局线程参数(Group)和本地线程参数(Local)。如图3所示，全局线程参数即全局存储空间，为整个磁场仿真加速计算空间的大小；本地线程参数为每个工作组的大小，即工作组矩阵的长×宽×高，本地存储空间为每个工作组的大小。

每个FDTD立方体网格具体位置坐标(I_x，I_y，I_z)通过下面三个式子得到：

I_x＝j_x×n_x+i_x         (17-1)

I_y＝j_y×n_y+i_y         (17-2)

I_z＝j_z×n_z+i_z         (17-3)

上面三式中j_x，j_y，j_z为边长与空间步长相等的立方体网格Δα³所在的工作组对应的空间位置；n_x，n_y，n_z为工作组沿三维的网格数；i_x，i_y，i_z为立方体网格Δα³在工作组内部所对应的具体位置坐标。

步骤4-4-3-3：定义总仿真计算空间位置坐标代号i000、i100、i010与i001，如下面式(18-1)到式(18-4)：

i000＝I_x×N_y×N_z/4+I_y×N_z/4+I_z        (18-1)

i100＝(I_x+1)×N_y×N_z/4+I_y×N_z/4+I_z    (18-2)

i010＝I_x×N_y×N_z/4+(I_y+1)×N_z/4+I_z    (18-3)

i001＝I_x×N_y×N_z/4+I_y×N_z/4+I_z+1      (18-4)

通过式(19-1)-(19-4)定义工作组空间位置坐标代号localindex、local100、local010和local001：

localindex＝i_x×n_y×n_z+i_y×n_z+i_z       (19-1)

local100＝(i_x+1)×n_y×n_z+i_y×n_z+i_z     (19-2)

local010＝i_x×n_y×n_z+(i_y+1)×n_z+i_z     (19-3)

local001＝i_x×n_y×n_z+i_y×n_z+i_z+1       (19-4)

利用步骤4-4-3-3的方法，可以充分利用float4型数据的并行计算特点，提高计算效率，并巧妙实现了FDTD算法相邻网格单元的电场相减。

步骤4-4-3-4：将GPU全局存储空间的电磁场数据传递到本地存储空间。

H_x-l[localindex]←H_x[i000]

H_y-l[localindex]←H_y[i000]

H_z-l[localindex]←H_z[i000]

E_x-l[localindex]←E_x[i000]

E_y-l[localindex]←E_y[i000]

E_z-l[localindex]←E_z[i000]

B_x-l[localindex]←B_x[i000]

B_y-l[localindex]←B_y[i000]

B_z-l[localindex]←B_z[i000]

其中，H_x、H_y、H_z分别表示全局存储空间的在x方向、y方向、z方向的磁场，H_x-l、H_y-l、H_z-l分别表示本地存储空间的在x方向、y方向、z方向的磁场，E_x、E_y、E_z分别表示全局存储空间的在x方向、y方向、z方向的电场，E_x-l、E_y-l、E_z-l分别表示本地存储空间的在x方向、y方向、z方向的电场，B_x、B_y、B_z分别表示全局存储空间的三维临时电磁参量，B_x-l、B_y-l、B_z-l分别表示本地存储空间的三维临时电磁参量。

步骤4-4-3-5：将三维磁场空间位置坐标的循环迭代转化为基于OpenCL的并行计算。具体设置方法如图4-图6所示。即将磁场仿真计算时，FDTD三维仿真计算空间的相邻电场之差转化为OpenCL所支持的本地存储空间内的float4型数据的相邻电场之差，在利用GPU仿真计算时，将充分利用GPU的流处理特性和float4型数据实现并行计算。

如图4所示，为y轴方向E_x和E_z数据转化过程。首先判断该网格处的坐标I_y是否小于FDTD的计算空间y轴坐标上限(N_y-1)。若否，则将E_x-010和E_z-010都设置为float4[1.0，1.0，1.0，1.0]。若是，则进一步判断本地存储空间的网格坐标i_y是否小于本地存储空间y轴纵坐标上限(n_y-1)，若否，则将E_x-010和E_z-010分别设置为E_x[i010]和E_z[i010]，若是，则将E_x-010和E_z-010分别设置为E_x-l[local010]和E_z-l[local010]，其中“E_x-010”和“E_z-010”的下角标010表示电场数据在y轴方向后移一位。

如图5所示，为z轴方向E_x和E_y数据转化过程。首先判断该网格处的纵坐标I_z是否小于FDTD的全局存储空间z轴坐标上限(N_z/4-1)；若否，则将E_x-001设置为float4(E_x-l[localindex].y，E_x-l[localindex].z，E_x-l[localindex].w，1.0)，将E_y-001设置为：float4(E_y-l[localindex].y，E_y-l[localindex].z，E_y-l[localindex].w，1.0)；若是，则进一步判断本地存储空间的网格坐标i_z是否小于本地存储空间z轴坐标上限(n_z-1)，若否，则将E_x-001和E_y-001分别设置为float4(E_x-l[localindex].y，E_x-l[localindex].z，E_x-l[localindex].w，E_x-l[i001].x)和float4(E_y-l[localindex].y，E_y-l[localindex].z，E_y-l[localindex].w，E_y-l[i001].x)，若是，则设置E_x-001为float4(E_x-l[localindex].y，E_x-l[localindex].z，E_x-l[localindex].w，E_x-l[local001].x)，设置E_y-001为float4(E_y-l[localindex].y，E_y-l[localindex].z，E_y-l[localindex].w，E_y-l[local001].x)。其中.x、.y、.z和.w依次代表float4型数据中的4个float型数据；E_x-001和E_y-001的下角标001表示电场数据在z轴方向后移一位。

如图6所示，为x轴方向E_y和E_z数据转化过程。首先判断该网格处的坐标I_x是否小于FDTD的全局存储空间x轴坐标上限(N_x-1)。若否，则将E_y-100和E_z-100都设置为float4[1.0，1.0，1.0，1.0]。若是，则进一步判断本地存储空间的网格坐标i_x是否小于本地存储空间x轴纵坐标上限(n_x-1)，若否，则将E_y-100和E_z-100分别设置为E_y[i100]和E_z[i100]；若是，则将E_y-100和E_z-100分别设置为E_y-l[local100]和E_z-l[local100]。E_y-100和E_z-100的下角标100表示电场数据在x轴方向后移一位。

步骤4-4-3-6：确定本地存储空间x方向的磁场H_x-l、y方向的磁场H_y-l、z方向的磁场H_z-l。

本地存储空间的x方向的临时电磁参量B_x-l如式(20)所示，得到的本地存储空间在x方向的磁场H_x-l如式(21)所示：

$B_{x-l}\left[\mathrm{localindex}\right]=\frac{2\mathrm{\ϵ}{\mathrm{\κ}}_y-{\mathrm{\σ}}_y\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}{B}_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\ϵ\Δt}}{2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}[\frac{{E_{y-001}}^{n-1}-E_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}}{\mathrm{\Δz}}-\frac{{E_{z-010}}^{n-1}-E_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}}{\mathrm{\Δy}}]$

$H_{x-l}\left[\mathrm{localindex}\right]=\left(\frac{2{\mathrm{\ϵ\κ}}_z-{\mathrm{\σ}}_z\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}\right)H_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\left[\frac{1}{(2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt})\mathrm{\μ}}\right]\left[(2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt})B_{x-l}\right[\mathrm{localindex}]-(2{\mathrm{\ϵ\κ}}_x-{\mathrm{\σ}}_x\mathrm{\Δt})B_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

本地存储空间的y方向的临时电磁参量B_y-l如式(22)所示，得到的本地存储空间在y方向的磁场H_y-l如式(23)所示：

$B_{y-l}\left[\mathrm{localindex}\right]=\frac{2{\mathrm{\ϵ\κ}}_z-{\mathrm{\σ}}_z\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}{B}_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\ϵ\Δt}}{2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}[\frac{{E_{z-100}}^{n-1}-E_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}}{\mathrm{\Δx}}-\frac{{E_{x-001}}^{n-1}-E_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}}{\mathrm{\Δz}}]$

$H_{y-l}\left[\mathrm{localindex}\right]=\left(\frac{2{\mathrm{\ϵ\κ}}_x-{\mathrm{\σ}}_x\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}\right)H_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\left[\frac{1}{({2\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt})\mathrm{\μ}}\right]\left[(2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt})B_{y-l}\right[\mathrm{localindex}]-(2{\mathrm{\ϵ\κ}}_y-{\mathrm{\σ}}_y\mathrm{\Δt})B_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

本地存储空间的z方向的临时电磁参量B_z-l如式(24)所示，得到的本地存储空间在z方向的磁场H_z-l如式(25)所示

$B_{z-l}\left[\mathrm{localindex}\right]=\frac{2{\mathrm{\ϵ\κ}}_x-{\mathrm{\σ}}_x\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}{B}_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\ϵ\Δt}}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}[\frac{{E_{x-010}}^{n-1}-E_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}}{\mathrm{\Δy}}-\frac{{E_{y-100}}^{n-1}-E_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}}{\mathrm{\Δx}}]$

$H_{z-l}\left[\mathrm{localindex}\right]=\left(\frac{2{\mathrm{\ϵ\κ}}_y-{\mathrm{\σ}}_y\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}\right)H_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\left[\frac{1}{(2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt})\mathrm{\μ}}\right]\left[(2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt})B_{z-l}\right[\mathrm{localindex}]-(2{\mathrm{\ϵ\κ}}_z-{\mathrm{\σ}}_z\mathrm{\Δt})B_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

步骤4-4-4：根据用户需要判断是否保存当前时间步下的磁场数据。如果计算的结果需要存储下来，则利用OpenCL库函数clEnqueueReadBuffer()将磁场H_x、H_y、H_z由GPU显存传递到CPU内存，以便于在内存中对磁场数据进行后处理，如果不需要则不进行此操作。

步骤4-5：执行基于OpenCL的电场计算内核函数。

步骤4-5-1：在时间步n＝1时，设置FDTD计算所用的内核参数，利用OpenCL标准库中clSetKernelArg()函数取得FDTD电场计算所用的三维电磁场数据变量以及三维临时电磁参量的存储地址。

步骤4-5-2：当利用GPU进行计算时，利用OpenCL标准库函数clEnqueueWriteBuffer()，将指令序列数据缓冲区中FDTD电场计算所用的的三维电磁场数据变量以及三维临时电磁参量传递到GPU显存。

步骤4-5-3：利用OpenCL标准库函数clEnqueueNDRangeKernel()激活GPU内核进行FDTD电场加速仿真。基于OpenCL内核函数的FDTD电场加速仿真如下步骤：

步骤4-5-3-1：定义全局线程参数(Group)和本地线程参数(Local)。全局参数为整个FDTD仿真空间的大小，即FDTD仿真空间总的长×宽×高；本地参数为每个工作组的大小，即工作组矩阵的长×宽×高。

每个FDTD网格具体位置坐标计算方式如下式(26-1)至式(26-3)所示。

I_x＝j_x×n_x+i_x    (26-1)

I_y＝j_y×n_y+i_y    (26-2)

I_z＝j_z×n_z+i_z    (26-3)

(26)式中j_x，j_y，j_z为该立方体网格Δα³所在的工作组对应的空间位置；n_x，n_y，n_z为工作组沿三维的网格数；i_x，i_y，i_z为立方体网格Δα³在工作组内部所对应的具体位置坐标。

步骤4-5-3-2：定义总的FDTD仿真空间位置坐标代号为式(27-1)至式(27-4)：

i000＝I_x×N_y×N_z/4+I_y×N_z/4+I_z        (27-1)

i100＝(I_x-1)×N_y×N_z/4+I_y×N_z/4+I_z    (27-2)

i010＝I_x×N_y×N_z/4+(I_y-1)×N_z/4+I_z    (27-3)

i001＝I_x×N_y×N_z/4+I_y×N_z/4+I_z-1      (27-4)

定义本地存储空间位置坐标代号为式(28-1)至式(28-4)

localindex＝i_x×n_y×n_z+i_y×n_z+i_z      (28-1)

local100＝(i_x-1)×n_y×n_z+i_y×n_z+i_z    (28-2)

local010＝i_x×n_y×n_z+(i_y-1)×n_z+i_z    (28-3)

local001＝i_x×n_y×n_z+i_y×n_z+i_z-1      (28-4)

利用4-5-3-2的方法，可以充分利用float4型数据的并行计算特点，提高计算效率，并巧妙实现了FDTD算法相邻网格单元的磁场相减。

步骤4-5-3-3：将GPU全局存储空间的电磁场数据传递到本地存储空间。

H_x-l[localindex]←H_x[i000]

H_y-l[localindex]←H_y[i000]

H_z-l[localindex]←H_z[i000]

E_x-l[localindex]←E_x[i000]

E_y-l[localindex]←E_y[i000]

E_z-l[localindex]←E_z[i000]

G_x-l[localindex]←G_x[i000]

G_y-l[localindex]←G_y[i000]

G_z-l[localindex]←G_z[i000]

F_x-l[localindex]←F_x[i000]

F_y-l[localindex]←F_y[i000]

F_z-l[localindex]←F_z[i000]

G_x、G_y、 G_z、F_x、F_y和F_z表示全局存储空间的三维临时电磁参量，G_x-l、G_y-l、G_z-l、F_x-l、F_y-l和F_z-l表示本地存储空间的三维临时电磁参量。全局存储空间的三维临时电磁参量如B_x，B_y，B_z，D_x，D_y，D_z，F_x，F_y，F_z，以及本地存储空间的三维临时电磁参量B_x-l，B_y-l，B_z-l，G_x-l，G_y-l，G_z-l，F_x-l，F_y-l和F_z-l都是在进行电磁场计算时的中间变量，没有物理意义，引入这些中间变量用于缩短公式长度。

步骤4-5-3-4：将三维电场空间位置坐标的循环迭代转化为基于OpenCL的并行计算。具体设置方法如图7-图9。即将电场计算时，FDTD三维计算空间的相邻磁场之差转化为OpenCL所支持的本地存储空间内的float4型数据的相邻磁场之差，在利用GPU计算时，将充分利用GPU的流处理特性和float4型数据实现并行计算。利用本方法的创新点是可以将传统FDTD的空间位置循环迭代转化为基于GPU流处理的并行计算。

如图7所示，为y轴方向H_x和H_z数据转化过程。首先判断该网格处的坐标I_y是否大于全局存储空间y轴坐标下限0。若否，则将H_x-010和H_z-010都设置为float4[1.0，1.0，1.0，1.0]。若是，则进一步判断本地存储空间的网格坐标i_y是否大于本地存储空间Y轴纵坐标下限0，若否，则将H_x-010和H_z-010分别设置为H_x[i010]和H_z[i010]，若是，则将H_x-010和H_z-010分别设置为H_x-l[local010]和H_z-l[local010]，其中，H_x-010和H_z-010的下角标010表示磁场数据在y轴方向前移一位。

如图8所示，z轴方向H_x和H_y数据转化过程。首先判断该网格处的纵坐标I_z是否大于FDTD的计算空间z轴坐标下限0。若否，设置H_x-001为float4(1.0f，H_x-l[localindex].x，H_x-l[localindex].y，H_x-l[localindex].z)，将H_y-001设置为：float4(1.0f，H_y-l[localindex].x，H_y-l[localindex].y，H_y-l[localindex].z)。若是，则进一步判断本地存储空间的网格坐标i_z是否大于本地存储空间Z轴坐标下限0。若否，则将H_x-001和H_y-001分别设置为float4(H_x-l[i001].w，H_x-l[localindex].x，H_x-l[localindex].y，H_x-l[localindex].z)和float4(H_y-l[i001].w，H_y-l[localindex].x，H_y-l[localindex].y，H_y-l[localindex].z)；若是，设置H_x-001为float4(H_x-l[local001].w，H_x-l[localindex].x，H_x-l[localindex].y，H_x-l[localindex].z)，设置H_y-001为float4(H_y-l[local001].w，H_y-l[localindex].x，H_y-l[localindex].y，H_y-l[localindex].z)。其中，H_x-001和H_y-001的下角标001表示磁场数据在z轴方向前移一位，.x、.y、.z和.w依次代表float4型数据中的4个float型数据。

如图9所示，x轴方向H_y和H_z数据转化过程。首先判断该网格处的坐标I_x是否大于FDTD的计算空间x轴坐标下限0。若否，则将H_y-100和H_z-100都设置为float4[1.0，1.0，1.0，1.0]。若是，则进一步判断本地存储空间的网格坐标i_x是否大于本地存储空间X轴纵坐标下限0，若否，则将H_y-100和H_z-100分别设置为H_y[i100]和H_z[i100]。若是，则将H_y-100和H_z-100分别设置为H_y-l[local100]和H_z-l[local100]。其中，H_y-100和H_z-100的下角标100表示磁场数据在x轴方向前移一位。

步骤4-5-3-5：计算本地存储空间的在x方向、y方向与z方向的电场。

本地存储空间的x方向的临时电磁参量F_x-l、G_x-l如式(29)与式(30)所示，得到的本地存储空间在x方向的电场E_x-l式(31)所示：

$F_{x-l}\left[\mathrm{localindex}\right]=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{F}_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}[\frac{H_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}-{H_{z-010}}^{n-1}}{\mathrm{\Δy}}-\frac{H_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}-{H_{y-001}}^{n-1}}{\mathrm{\Δz}}]---\left(29\right)$

$G_{x-l}\left[\mathrm{localindex}\right]=\frac{2\mathrm{\ϵ}{\mathrm{\κ}}_y-{\mathrm{\σ}}_y\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}{G}_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\ϵ}}{{2\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}\left[F_{x-l}\right[\mathrm{localindex}]-F_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

$E_{x-l}\left[\mathrm{localindex}\right]=\frac{{2\mathrm{\ϵ\κ}}_z-{\mathrm{\σ}}_z\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}{E}_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{1}{2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}\left[(2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt})G_{x-l}\right[\mathrm{localindex}]-({2\mathrm{\ϵ\κ}}_x-{\mathrm{\σ}}_x\mathrm{\Δt})G_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

本地存储空间的y方向的临时电磁参量F_y-l、G_y-l如式(32)与式(33)所示，得到的本地存储空间在y方向的电场E_y-l如式(34)所示：

$F_{y-l}\left[\mathrm{localindex}\right]=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{F}_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}[\frac{H_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}-{H_{x-001}}^{n-1}}{\mathrm{\Δz}}-\frac{H_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}-{H_{z-100}}^{n-1}}{\mathrm{\Δx}}]---\left(32\right)$

$G_{y-l}\left[\mathrm{localindex}\right]=\frac{2{\mathrm{\ϵ\κ}}_z-{\mathrm{\σ}}_z\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}{G}_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\ϵ}}{{2\mathrm{\ϵ\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt}}\left[F_{y-l}\right[\mathrm{localindex}]-F_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

$E_{y-l}\left[\mathrm{localindex}\right]=\frac{2{\mathrm{\ϵ\κ}}_x-{\mathrm{\σ}}_x\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}{E}_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{1}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}\left[(2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt})G_{y-l}\right[\mathrm{localindex}]-({2\mathrm{\ϵ\κ}}_y-{\mathrm{\σ}}_y\mathrm{\Δt})G_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

本地存储空间的z方向的临时电磁参量F_z-l、G_z-l如式(35)与式(36)所示，得到的本地存储空间在z方向的电场E_Z-l如式(37)所示：

$F_{z-l}\left[\mathrm{localindex}\right]=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{F}_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}[\frac{H_{y-l}{\left[\mathrm{localindex}\right]}^{n-1}-{H_{y-100}}^{n-1}}{\mathrm{\Δx}}-\frac{H_{x-l}{\left[\mathrm{localindex}\right]}^{n-1}-{H_{x-010}}^{n-1}}{\mathrm{\Δy}}]---\left(35\right)$

$G_{z-l}\left[\mathrm{localindex}\right]=\frac{2{\mathrm{\ϵ\κ}}_x-{\mathrm{\σ}}_x\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}{G}_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{2\mathrm{\ϵ}}{2{\mathrm{\ϵ\κ}}_x+{\mathrm{\σ}}_x\mathrm{\Δt}}\left[F_{z-l}\right[\mathrm{localindex}]-F_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

$E_{z-l}\left[\mathrm{localindex}\right]=\frac{2{\mathrm{\ϵ\κ}}_y-{\mathrm{\σ}}_y\mathrm{\Δt}}{2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}{E}_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}$

$+\frac{1}{2{\mathrm{\ϵ\κ}}_y+{\mathrm{\σ}}_y\mathrm{\Δt}}\left[(2\mathrm{\ϵ}{\mathrm{\κ}}_z+{\mathrm{\σ}}_z\mathrm{\Δt})G_{z-l}\right[\mathrm{localindex}]-(2{\mathrm{\ϵ\κ}}_z-{\mathrm{\σ}}_z\mathrm{\Δt})G_{z-l}{\left[\mathrm{localindex}\right]}^{n-1}]$

步骤4-5-4：根据是否需要对该电场数据进行后处理来判断是否需要保存当前时间步下的电场数据。如果需要保存，则利用OpenCL标准库函数clEnqueueReadBuffer()将电场E_x、E_y、E_z由GPU显存传递到CPU内存，以便于在内存中对电场数据进行后处理。

步骤4-6：保存所需要的电磁场数据。

步骤4-6-1：保存时间步n₀下的空间电磁场数据。

以平面的形式存储电磁场，设置所需要存储的计算空间内的平面位置。平面电磁场的存储方式可以分为将E_x，E_y，E_z，H_x，H_y和H_z按垂直于x-y面、垂直于y-z面和垂直于x-z面分别储存，即存储18个电磁场数据文件，如图10所示。然后将CPU内存中的电磁场数据存储到文件中，以便分析使用。

步骤4-6-2：保存所关心的位置点的电磁场在时间步n₁～n₂情况下的电场与磁场数据随时间的变化情况，以便分析使用。

步骤5：释放OpenCL存储资源。当FDTD计算和数据存储完成后，利用OpenCL标准库函数clReleaseMemObject()释放显存中存储的数据。

步骤6：仿真完成后，释放内存资源。当GPU加速FDTD完成后，释放CPU内存中的数据。

为了验证本发明方法加速FDTD电磁场方法效率，进行了如下的数值仿真试验。GPU型号采用AMD HD4850，CPU采用同等时期的双核Intel Dual E2180。仿真系统为Windows XP 32-bit操作系统。

如表1与图11所示，在固定迭代时间步数n为1000，改变FDTD的仿真空间条件下GPU和CPU仿真时间的对比，从图11中(a)可看出，基于GPU的仿真时间明显低于基于CPU的仿真时间，随着FDTD仿真空间的增加，基于CPU仿真时间越来越长，基于GPU仿真的时间并没有基于CPU仿真的时间增长的比例快，从图11中(b)可以看出，基于CPU仿真的时间是基于GPU仿真时间的5倍以上。因此，可见基于OpenCL的GPU加速三维电磁场FDTD仿真效果明显。

表1.固定迭代时间步数1000，改变FDTD仿真空间条件下的GPU和CPU的仿真时间表

注：小数点后第四位四舍五入

如表2与图12所示，在固定FDTD仿真空间尺寸为120×120×120网格，改变迭代时间步数条件下GPU和CPU仿真时间的对比，从图12中的(a)可以看出基于GPU仿真的时间明显低于基于CPU仿真的时间，从图12中(b)可以看出，基于CPU仿真的时间是基于GPU仿真时间的5倍以上。

表2.固定FDTD仿真空间尺寸120×120×120网格数，改变迭代时间步数条件下GPU和CPU的仿真时间表

注：小数点后第四位四舍五入

如表3与图13所示，在固定迭代时间步数1000，改变FDTD仿真空间大小的条件下基于OpenCL的双核并行CPU仿真时间与普通CPU仿真时间对比，如图13中(a)与(b)可以看出基于OpenCL的双核并行CPU的仿真时间明显低于基于普通CPU的仿真时间。

表3.固定迭代时间步数1000，改变FDTD仿真空间条件下基于OpenCL的双核并行CPU与普通CPU的仿真时间

注：小数点后第四位四舍五入

Claims (8)

1.一种基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，该方法具体包括以下步骤：

步骤1：设置FDTD仿真参数，为所述仿真参数动态申请内存空间；

步骤2：构建三维的电磁模型，对该电磁模型进行网格剖分；

步骤3：初始化OpenCL的计算参数，包括：获得FDTD仿真空间的大小、设置仿真所用处理器类型、创建上下文、获得工作组的维数及大小信息、创建指令序列及数据缓存区、装载FDTD加速程序文件、声明磁场仿真计算的内核函数以及电场仿真计算的内核函数；

步骤4：基于OpenCL对三维电磁模型进行FDTD加速仿真；

步骤5：释放显存中存储的数据；

步骤6：仿真完成后，释放CPU内存资源；

其中，OpenCL表示开放计算语言，GPU表示图形处理器，FDTD表示时域有限差分。

2.根据权利要求1所述的基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，步骤1具体包括：

步骤1-1：设置电磁场仿真参数，包括如下步骤：

步骤1-1-1：设置空间步长Δα：

$\mathrm{\Δ\α}=\frac{{\mathrm{\λ}}_m}{10},$ α＝x，y，z    (1)

λ_m为激励源的最高频率所对应的波长；

步骤1-1-2：设置时间步长Δt，时间步长Δt的设置满足式(2)：

$\mathrm{\Δt}=\frac{0.9}{c\sqrt{\frac{1}{\mathrm{\Δ}{x}^2}+\frac{1}{\mathrm{\Δ}{y}^2}+\frac{1}{\mathrm{\Δ}{z}^2}}}---\left(2\right)$

c为真空中的光速；

步骤1-1-3：设置UPML吸收边界：设置UPML吸收边界在x-y-z三维方向的网格层数n_x-PML1，n_x-PML2，n_y-PML1，n_y-PML2，n_z-PML1和n_z-PML2，以及吸收边界参数σ和κ；

步骤1-1-4：激励源设置，包括电流源设置与偶极子源设置；

电流源设置：在麦克斯韦磁场旋度方程的右侧添加电流密度项，如式(3)所示，根据实际情况，选择在x-y-z三维方向的其中一个方向上施加该电流源：

$\▿\×\stackrel{\→}{H}=\frac{\∂\stackrel{\→}{D}}{\∂t}+\stackrel{\→}{J}---\left(3\right)$

其中，

为电流密度矢量，为磁场矢量，

为电位移矢量，t为时间；

偶极子源设置：在麦克斯韦磁场旋度方程的右侧施加电偶极子项，如式(4)所示，根据实际情况，选则在x-y-z三维方向的其中一个方向上施加该偶极子源：

$\▿\×\stackrel{\→}{H}=\mathrm{\ϵ}\frac{\∂\stackrel{\→}{E}}{\∂t}+\frac{1}{{\mathrm{\Δ\α}}^3}\frac{e\stackrel{\→}{p}}{\mathrm{dt}}---\left(4\right)$

其中，

表示偶极子源矢量，

表示电场矢量，ε表示介电常数；

步骤1-2：动态申请内存空间：采用动态分配方式，连续存储数据，对于三维电磁场数据变量以及由UPML引入的三维临时电磁参量，其寻址方式为：

                 φ[i×N_yz+j×N_z+k]，N_yz＝N_y×N_z    (5)

其中i，j和k为分别代表三维坐标系x轴、y轴和z轴上所对应的空间位置点，N_y和N_z为y轴和z轴上所对应的网格数。

3.根据权利要求1所述的基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，所述的步骤2具体包括以下两步：

步骤2-1：简单的电磁模型的构建及网格剖分：直接对该电磁模型建立三维几何模型，并进行网格剖分；

步骤2-2：复杂的电磁模型的构建及网格自动剖分：利用建模软件Auto CAD构建三维几何模型，或读入已经建立的三维几何模型，不同电磁材料部件用图层来区分，然后利用二次开发工具Object ARX开发自动剖分程序输出有效网格位置及其对应电磁材料参数，生成几何-电磁参数描述文件，然后调用该描述文件，将剖分网格和电磁参数读入计算内存。

4.根据权利要求1所述的基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，所述的步骤3具体为：

步骤3-1：确定基于GPU的FDTD仿真空间的大小：N_x×N_y×N_z/4；

步骤3-2：设置仿真所用处理器类型：从OpenCL标准库中选取参数CL_DEVICE_TYPE_GPU进行设置；所述的CL_DEVICE_TYPE_GPU表示仿真所用处理器为GPU；

步骤3-3：利用OpenCL标准库函数clCreateContextFromType()创建上下文，并利用OpenCL标准库函数clGetContextInfo()获得上下文信息；

步骤3-4：利用OpenCL标准库函数clGetDeviceInfo()获得OpenCL的工作组的维数及大小信息；

步骤3-5：利用OpenCL标准库函数clCreateCommandQueue()创建指令序列，并创建三维电磁场数据变量以及由UPML引入的三维临时电磁参量的数据缓冲区；

步骤3-6：利用OpenCL标准库函数clCreateProgramWithSource()装载基于OpenCL的FDTD加速程序文件，并利用OpenCL标准库函数clBuildProgram()将加速程序文件转化为GPU内核可执行的程序文件；

步骤3-7：创建声明基于OpenCL的磁场仿真计算的内核函数；

步骤3-8：创建声明基于OpenCL的电场仿真计算的内核函数。

5.根据权利要求1或2所述的基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，所述的步骤4具体包括以下步骤：

步骤4-1：设置时间迭代步数n，则总的迭代时间为t＝n×Δt；

步骤4-2：为每一个迭代步设置其对应的开始时间：t₀+Δt×(n-1)，其中t₀为仿真起始时刻；

步骤4-3：按照步骤1中的激励源设置方法设置激励源；

步骤4-4：执行基于OpenCL的FDTD磁场仿真计算内核函数：

步骤4-4-1：在时间步n＝1时，利用OpenCL标准库函数clSetKernelArg()取得FDTD磁场仿真计算所用的三维电磁场数据变量以及三维临时电磁参量的存储地址；

步骤4-4-2：当基于GPU进行仿真时，利用OpenCL标准库函数clEnqueueWriteBuffer()，将指令序列数据缓冲区中FDTD磁场仿真计算所用的三维电磁场数据变量以及三维临时电磁参量传递到GPU显存；

步骤4-4-3：利用OpenCL标准库函数clEnqueueNDRangeKernel()激活GPU内核，进行FDTD磁场加速仿真；

步骤4-4-4：如果需要对仿真后的磁场数据进行存储，则利用OpenCL标准库函数clEnqueueReadBuffer()将x方向、y方向、z方向的磁场H_x、H_y、H_z由GPU显存传递到CPU内存，以便在内存中对磁场数据进行后处理；

步骤4-5：执行基于OpenCL的FDTD电场仿真计算内核函数：

步骤4-5-1：在时间步n＝1时，利用OpenCL标准库函数clSetKernelArg()取得FDTD电场仿真计算所用的三维电磁场数据变量以及三维临时电磁参量的存储地址；

步骤4-5-2：当基于GPU进行仿真时，利用OpenCL标准库函数clEnqueueWriteBuffer()，将指令序列数据缓冲区中FDTD电场计算所用的三维电磁场数据变量以及三维临时电磁参量传递到GPU显存；

步骤4-5-3：利用OpenCL标准库函数clEnqueueNDRangeKernel()激活GPU内核进行FDTD电场加速仿真；

步骤4-5-4：如果需要对仿真后的电场数据进行存储，则利用OpenCL标准库函数clEnqueueReadBuffer()将电场E_x、E_y、E_z由GPU显存传递到CPU内存，以便于在内存中对电场数据进行后处理；

步骤4-6：保存所需要的电磁场数据：保存时间步n₀下的空间电磁场数据，以及保存所关心的位置点的电磁场在时间步n₁～n₂下的电磁场随时间变化情况的数据。

6.根据权利要求5所述的基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，步骤4-4-3所述的磁场仿真加速计算的方法如下：

步骤4-4-3-1：设置工作组大小为：n_x×n_y×n_z＝2×2×30；

步骤4-4-3-2：设置全局存储空间为整个FDTD仿真空间的大小，本地存储空间为每个工作组的大小，每个FDTD立方体网格具体位置坐标(I_x，I_y，I_z)为：

$\left\{\begin{array}{c}{I}_x=j_x\×n_x+i_x\\ I_y=j_y\×n_y+i_y\\ I_z=j_z\×n_z+i_z\end{array}\right.---\left(6\right)$

其中，j_x，j_y，j_z为边长与空间步长相等的立方体网格Δα³所在的工作组对应的空间位置；n_x，n_y，n_z为工作组沿三维的网格数；i_x，i_y，i_z为立方体网格Δα³在工作组内部所对应的具体位置坐标；

步骤4-4-3-3：通过式(7)确定全局存储空间的位置坐标代号i000、i100、i010与i001：

$\left\{\begin{array}{c}i000=I_x\×N_y\×N_z/4+I_y\×N_z/4+I_z\\ i100=(I_x+1)\×N_y\×N_z/4+I_y\×N_z/4+I_z\\ i010=I_x\×N_y\×N_z/4+(I_y+1)\×N_z/4+I_z\\ i001=I_x\×N_y\×N_z/4+I_y\×N_z/4+I_z+1\end{array}\right.---\left(7\right)$

N_y和N_z为整个FDTD仿真空间y轴和z轴方向所对应的网格；

通过式(8)确定本地存储空间的位置坐标代号localindex、local100、local010和local001：

$\left\{\begin{array}{c}\mathrm{localindex}=i_x\×n_y\×n_z+i_y\×n_z+i_z\\ \mathrm{local}100=(i_x+1)\×n_y\×n_z+i_y\×n_z+i_z\\ \mathrm{local}010=i_x\×n_y\×n_z+(i_y+1)\×n_z+i_z\\ \mathrm{local}001=i_x\×n_y\×n_z+i_y\×n_z+i_z+1\end{array}\right.---\left(8\right)$

步骤4-4-3-4：将GPU全局存储空间的电磁场数据及相关的三维临时电磁参量传递到本地存储空间；

步骤4-4-3-5：将三维磁场空间位置坐标的循环迭代转化为基于OpenCL的并行计算；

(1)y轴方向E_x和E_z的数据转化过程：首先判断该网格处的坐标I_y是否小于全局存储空间Y轴坐标上限(N_y-1)，若否，则将E_x-010和E_z-010都设置为float4[1.0，1.0，1.0，1.0]；若是，则进一步判断本地存储空间的网格坐标i_y是否小于本地存储空间y轴纵坐标上限(n_y-1)，若否，则将E_x-010和E_z-010分别设置为E_x[i010]和E_z[i010]，若是，则将E_x-010和E_z-010分别设置为E_x-l[local010]和E_z-l[local010]；其中E_x-010和E_z-010的下角标010表示电场数据在y轴方向后移一位；

(2)z轴方向E_x和E_y，的数据转化过程：首先，判断该网格处的纵坐标I_z是否小于全局存储空间z轴坐标上限(N_z/4-1)，若不是，则将E_x-001设置为float4(E_x-l[localindex].y，E_x-l[localindex].z，E_x-l[localindex].w，1.0)，将E_y-001设置为float4(E_y-l[localindex].y，E_y-l[localindex].z，E_y-l[localindex].w，1.0)；若是，则进一步判断本地存储空间的网格坐标i_z是否小于本地存储空间z轴坐标上限(n_z-1)，若否，则将E_x-001和E_y-001分别设置为float4(E_x-l[localindex].y，E_x-l[localindex].z，E_x-l[localindex].w，E_x-l[i001].x)和float4(E_y-l[localindex].y，E_y-l[localindex].z，E_y-l[localindex].w，E_y-l[i001].x)，若是，则设置E_x-001为float4(E_x-l[localindex].y，E_x-l[localindex].z，E_x-l[localindex].w，E_x-l[local001].x)，设置E_y-001为float4(E_y-l[localindex].y，E_y-l[localindex].z，E_y-l[localindex].w，E_y-l[local001].x)；其中.x、.y、.z和.w依次表示float4型数据中4个float型数据，E_x-001和E_y-001的下角标001表示电场数据在z轴方向后移一位；

(3)x轴方向E_y和E_z数据转化过程：首先判断该网格处的坐标I_x是否小于全局存储空间x轴坐标上限(N_x-1)，若否，则将E_y-100和E_z-100都设置为float4[1.0，1.0，1.0，1.0；若是，则进一步判断本地存储空间的网格坐标i_x是否小于本地存储空间x轴纵坐标上限(n_x-1)，若否，则将E_y-100和E_z-100分别设置为E_y[i100]和E_z[i100]，若是，则将E_y-100和E_z-100分别设置为E_y-l[local100]和E_z-l[local100]；其中E_y-100和E_z-100的下角标100表示电场数据在x轴方向后移一位；

步骤4-4-3-6：确定x方向的磁场H_x、y方向的磁场H_y与z方向的磁场H_z。

7.根据权利要求5所述的基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，步骤4-5-3所述的电场加速计算方法如下：

步骤4-5-3-1：设置全局存储空间为整个FDTD仿真空间的大小，本地存储空间为每个工作组的大小，每个FDTD立方体网格具体位置坐标(I_x，I_y，I_z)为：

$\left\{\begin{array}{c}{I}_x=j_x\×n_x+i_x\\ I_y=j_y\×n_y+i_y\\ I_z=j_z\×n_z+i_z\end{array}\right.---\left(9\right)$

其中，j_x，j_y，j_z为边长与空间步长相等的立方体网格Δα³所在的工作组对应的空间位置；n_x，n_y，n_z为工作组沿三维的网格数；i_x，i_y，i_z为立方体网格Δα³在工作组内部所对应的具体位置坐标；

步骤4-5-3-2：通过式(10)确定全局存储空间的位置坐标代号i000、i100、i010与i001：

$\left\{\begin{array}{c}i000=I_x\×N_y\×N_z/4+I_y\×N_z/4+I_z\\ i100=(I_x+1)\×N_y\×N_z/4+I_y\×N_z/4+I_z\\ i010=I_x\×N_y\×N_z/4+(I_y+1)\×N_z/4+I_z\\ i001=I_x\×N_y\×N_z/4+I_y\×N_z/4+I_z+1\end{array}\right.---\left(10\right)$

N_y和N_z为整个FDTD仿真空间Y轴和Z轴方向所对应的网格；

通过式(11)确定本地存储空间的位置坐标代号localindex、local100、local010和local001：

$\left\{\begin{array}{c}\mathrm{localindex}=i_x\×n_y\×n_z+i_y\×n_z+i_z\\ \mathrm{local}100=(i_x+1)\×n_y\×n_z+i_y\×n_z+i_z\\ \mathrm{local}010=i_x\×n_y\×n_z+(i_y+1)\×n_z+i_z\\ \mathrm{local}001=i_x\×n_y\×n_z+i_y\×n_z+i_z+1\end{array}\right.---\left(11\right)$

步骤4-5-3-3：将GPU全局存储空间的电磁场数据及相关的三维临时电磁参量传递到本地存储空间；

步骤4-5-3-4：将三维电场空间位置坐标的循环迭代转化为基于OpenCL的并行计算；

(1)y轴方向H_x和H_z数据转化过程：首先判断该网格处的坐标I_y是否大于全局存储空间Y轴坐标下限0，若否，则将H_x-010和H_z-010都设置为float4[1.0，1.0，1.0，1.0]；若是，则进一步判断本地存储空间的网格坐标i_y是否大于本地存储空间y轴纵坐标下限0，若否，则将H_x-010和H_z-010分别设置为H_x[i010]和H_z[i010]，若是，则将H_x-010和H_z-010分别设置为H_x-l[local010]和H_z-l[local010]；

(2)z轴方向H_x和H_y数据转化过程：首先判断该网格处的纵坐标I_z是否大于全局存储空间z轴坐标下限0，若否，则将H_x-001设置为：float4(1.0f，H_x-l[localindex].x，H_x-l[localindex].y，H_x-l[localindex].z)，将H_y-001设置为：float4(1.0f，H_y-l[localindex].x，H_y-l[localindex].y，H_y-l[localindex].z)；若是，进一步判断本地存储空间的网格坐标i_z是否大于本地存储空间z轴坐标下限0，若否，则将H_x-001设置为float4(H_x-l[i001].w，H_x-l[localindex].x，H_x-l[localindex].y，H_x-l[localindex].z)，将H_y-001设置为float4(H_y-l[i001].w，H_y-l[localindex].x，H_y-l[localindex].y，H_y-l[localindex].z)，若是，设置H_x-001为float4(H_x-l[local001].w，H_x-l[localindex].x，H_x-l[localindex].y，H_x-l[localindex].z)，置H_y-001为float4(H_y-l[local001].w，H_y-l[localindex].x，H_y-l[localindex].y，H_y-l[localindex].z)；其中H_x-001和H_y-001的下角标001表示磁场数据在z轴方向前移一位，.x、.y、.z和.w依次表示float4型数据中4个float型数据；

(3)x轴方向H_y和H_z数据转化过程：首先判断该网格处的坐标I_x是否大于全局存储空间x轴坐标下限0，若否，则将H_y-100和H_z-100都设置为float4[1.0，1.0，1.0，1.0]；若是，则进一步判断本地存储空间的网格坐标i_x是否大于本地存储空间x轴纵坐标下限0，若否，则将H_y-100和H_z-100分别设置为H_y[i100]和H_z[i100]，若是，则将H_y-100和H_z-100分别设置为H_y-l[local100]和H_z-l[local100]；其中，H_y-100和H_z-100的下角标100表示磁场数据在x轴方向前移一位；

步骤4-5-3-5：确定x方向的电场E_x、y方向的电场E_y与z方向的电场E_z。

8.根据权利要求1所述的一种基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，所述的步骤5具体是利用OpenCL标准库函数clReleaseMemObject()来释放显卡中存储的数据。

Images