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

摘要

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

主权项

1.一种基于OpenCL的GPU加速三维时域有限差分电磁场仿真的方法，其特征在于，该方法具体包括以下步骤：步骤1：设置FDTD仿真参数，为所述仿真参数动态申请内存空间；步骤1具体包括：步骤1-1：设置电磁场仿真参数，包括如下步骤：步骤1-1-1：设置空间步长Δα：$\mathrm{\Delta \alpha }=\frac{{\lambda }_m}{10},\alpha =x,y,z---\left(1\right)$λ_m为激励源的最高频率所对应的波长;步骤1-1-2：设置时间步长Δt，时间步长Δt的设置满足式(2)：$\mathrm{\Delta t}=\frac{0.9}{c\sqrt{\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}+\frac{1}{\Delta 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三维方向的其中一个方向上施加该电流源：$▿\times \overrightarrow{H}=\frac{\partial \overrightarrow{D}}{\partial t}+\overrightarrow{J}---\left(3\right)$其中，▽表示旋度符号，为电流密度矢量，为磁场矢量，为电位移矢量，t为时间；偶极子源设置：在麦克斯韦磁场旋度方程的右侧施加电偶极子项，如式(4)所示，根据实际情况，选择在x-y-z三维方向的其中一个方向上施加该偶极子源：$▿\times \overrightarrow{H}=ϵ\frac{\partial \overrightarrow{E}}{\partial t}+\frac{1}{\Delta {\alpha }^3}\frac{d\overrightarrow{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轴上所对应的网格数；步骤2：构建三维的电磁模型，对该电磁模型进行网格剖分；步骤3：初始化OpenCL的计算参数，包括：获得FDTD仿真空间的大小、设置仿真所用处理器类型、创建上下文、获得工作组的维数及大小信息、创建指令序列及数据缓存区、装载FDTD加速程序文件、声明磁场仿真计算的内核函数以及电场仿真计算的内核函数；步骤4：基于OpenCL对三维电磁模型进行FDTD加速仿真；步骤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-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\times n_x+i_x\\ I_y=j_y\times n_y+i_y\\ I_z=j_z\times 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\times N_y\times N_z/4+I_y\times N_z/4+I_z\\ i100=(I_x+1)\times N_y\times N_z/4+I_y\times N_z/4+I_z\\ i010=I_x\times N_y\times N_z/4+(I_y+1)\times N_z/4+I_z\\ i001=I_x\times N_y\times N_z/4+I_y\times 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\times n_y\times n_z+i_y\times n_z+i_z\\ \mathrm{local}100=(i_x+1)\times n_y\times n_z+i_y\times n_z+i_z\\ \mathrm{local}010=i_x\times n_y\times n_z+(i_y+1)\times n_z+i_z\\ \mathrm{local}001=i_x\times n_y\times n_z+i_y\times 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；步骤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-3包括如下步骤：步骤4-5-3-1：设置全局存储空间为整个FDTD仿真空间的大小，本地存储空间为每个工作组的大小，每个FDTD立方体网格具体位置坐标(I_x,I_y，I_z)为：$\left\{\begin{array}{c}{I}_x=j_x\times n_x+i_x\\ I_y=j_y\times n_y+i_y\\ I_z=j_z\times 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\times N_y\times N_z/4+I_y\times N_z/4+I_z\\ i100=(I_x-1)\times N_y\times N_z/4+I_y\times N_z/4+I_z\\ i010=I_x\times N_y\times N_z/4+(I_y-1)\times N_z/4+I_z\\ i001=I_x\times N_y\times N_z/4+I_y\times 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\times n_y\times n_z+i_y\times n_z+i_z\\ \mathrm{local}100=(i_x-1)\times n_y\times n_z+i_y\times n_z+i_z\\ \mathrm{local}010=i_x\times n_y\times n_z+(i_y-1)\times n_z+i_z\\ \mathrm{local}001=i_x\times n_y\times n_z+i_y\times 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]；其中，H_x-010和H_z-010的下角标010表示磁场数据在y轴方向前移一位；（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；步骤4-5-4：如果需要对仿真后的电场数据进行存储，则利用OpenCL标准库函数clEnqueueReadBuffer()将电场E_x、E_y、E_z由GPU显存传递到CPU内存，以便于在内存中对电场数据进行后处理；步骤4-6：保存所需要的电磁场数据:保存时间步n₀下的空间电磁场数据，以及保存所关心的位置点的电磁场在时间步n₁～n₂下的电磁场随时间变化情况的数据；步骤5：释放显存中存储的数据；步骤6：仿真完成后，释放CPU内存资源；其中，OpenCL表示开放计算语言，GPU表示图形处理器，FDTD表示时域有限差分。