[发明专利]一种基于GPU并行的LeapfrogADI-FDTD方法

摘要

本发明属于计算电磁学与高性能计算领域，具体为一种基于GPU并行的Leapfrog ADI‑FDTD方法。本发明采用Yee网格，通过将GPU并行计算技术应用到含有CPML吸收边界的Leapfrog ADI‑FDTD算法中，形成了基于GPU并行加速的含有CPML吸收边界的Leapfrog ADI‑FDTD算法，在保持Leapfrog ADI‑FDTD算法与CPML吸收边界的优良特性的同时，极大地提高了计算效率，大大节约了计算时间，减少了数值模拟实验的时间周期和成本，达到了相关科学研究或器件设计的高效需求。

主权项

一种基于GPU并行的Leapfrog ADI‑FDTD方法，具体步骤如下：步骤1、将求解区域使用Yee网格划分；对CPU主机端与GPU上的数据进行初始化，即在CPU主机端定义结构体数组，大小为模拟求解区域的总网格数，用于存放GPU返回的各个网格点的电磁场计算结果，并将各个网格的电磁场场值赋值为0；同时，在GPU上的全局内存中定义同样大小的结构体数组，用于存储电磁场计算结果，并将各个网格的电磁场场值赋值为0；步骤2、数据预处理，在GPU全局内存上开辟空间，并利用下述推导的计算式，计算得到矩阵方程的系数矩阵；含有CPML的Leapfrog ADI‑FDTD算法计算的场量包括电场的三个分量Ex、Ey、Ez，磁场的三个分量Hx、Hy、Hz；其中，Ex的计算公式如式(1)、(2)、(3)所示，Hx的计算公式如式(4)、(5)、(6)所示：[1-Δt24μϵδykyj-1δy]Exn+12(i+12,j,k)=[1-Δt24μϵδykyj-1δy]Exn-12(i+12,j,k)+Δtϵ[1kyjδyHzn(i+12,j,k)-1kzkδzHyn(i+12,j,k)+ψexyn(i+12,j,k)-ψexzn(i+12,j,k)]---(1)]]>ψexyn(i+12,j,k)=byψexyn-1(i+12,j,k)+ay(Hzn(i+12,j+12,k)-Hzn(i+12,j-12,k)Δy)---(2)]]>ψexzn(i+12,j,k)=bzψexzn-1(i+12,j,k)+az(Hyn(i+12,j,k+12)-Hyn(i+12,j,k-12)Δz)---(3)]]>[1-Δt24μϵδykyj-1δy]Hxn+1(i,j+12,k+12)=[1-Δt24μϵδykyj-1δy]Hxn(i,j+12,k+12)+Δtμ[1kzkδzEyn+12(i,j+12,k+12)-1kyjδyEzn+12(i,j+12,k+12)+ψhxzn+12(i,j+12,k+12)-ψhxyn+12(i,j+12,k+12)]---(4)]]>ψhxzn+12(i,j+12,k+12)=bzψhxzn-12(i,j+12,k+12)+az(Eyn+12(i,j+12,k+1)-Eyn+12(i,j+12,k)Δz)---(5)]]>ψhxyn+12(i,j+12,k+12)=byψhxyn-12(i,j+12,k+12)+ay(Ezn+12(i,j+1,k+12)-Ezn+12(i,j,k+12)Δy)---(6)]]>其中，μ为磁导率，ε为介电常数，ks为坐标伸缩系数s＝x,y,z，Δt为单位时间步长，Δx、Δy、Δz分别为沿x、y、z方向的单位网格长度，as与bs为中间计算结果，二者的计算是为了得到ψ，其计算式如式(21)(20)所示，ψ为方程右端常数项的中间结果，i、j、k为模拟区域中网格在x、y、z方向上的编号，n为当前时间步长；其他方向的电场Ey、Ez与磁场Hy、Hz分量计算式同理可得；由于在Yee网格中，存在半网格位置，将地址映射到整数网格点，将(1)～(6)式进行地址映射与化简后我们得到Ex的计算式如式(7)～(12)所示，Hx的计算式如式(13)～(18)所示：AExExn+12(i,j-1,k)+BExExn+12(i,j,k)+CExExn+12(i,j+1,k)=AExExn-12(i,j-1,k)+BExExn-12(i,j,k)+CExExn-12(i,j+1,k)+Δt[(Hzn(i,j,k)-Hzn(i,j-1,k))/Δykyj-(Hyn(i,j,k)-Hyn(i,j,k-1))/Δzkzk+ψexyn(i,j,k)-ψexzn(i,j,k)](7)]]>其中：AEx=-(Δt2)2μz-1kyj-1(1Δy)2---(8)]]>BEx=(ϵx+2(Δt2)2μz-1kyj-1(1Δy)2)---(9)]]>CEx＝AEx  (10)ψexyn(i,j,k)=byψexyn-1(i,j,k)+ay(Hzn(i,j,k)-Hzn(i,j-1,k)Δy)---(11)]]>ψexzn(i,j,k)=bzψexzn-1(i,j,k)+az(Hyn(i,j,k)-Hyn(i,j,k-1)Δz)---(12)]]>AHxHxn+1(i,j-1,k)+BHxHxn+1(i,j,k)+CHxHxn+1(i,j+1,k)=AHxHxn(i,j-1,k)+BHxHxn(i,j,k)+CHxHxn(i,j+1,k)+Δt[(Eyn+12(i,j,k+1)-Eyn+12(i,j,k))/Δzkzk-(Ezn+12(i,j+1,k)-Ezn+12(i,j,k))/Δykyj+ψhxzn+12(i,j,k)-ψhxyn+12(i,j,k)]---(13)]]>其中：AHx=-(Δt2)2ϵz-1kyj-1(1Δy)2---(14)]]>BHx=(μx+2(Δt2)2ϵz-1kyj-1(1Δy)2)---(15)]]>CHx＝AHx  (16)ψhxzn+12(i,j,k)=bzψhxzn-12(i,j,k)+az(Eyn+12(i,j,k+1)-Eyn+12(i,j,k)Δz)---(17)]]>ψhxyn+12(i,j,k)=byψhxyn-12(i,j,k)+ay(Ezn+12(i,j+1,k)-Ezn+12(i,j,k)Δy)---(18)]]>其中，AEx、BEx、CEx为Ex计算方程中左端的系数,AHx、BHx、CHx为Hx计算方程中左端的系数；其他方向电场分量与磁场分量计算式同理可得；其中，方程中特定的变量求解式如下：(s＝x,y,z)  (19)bs=e-(σsks+αs)(Δtϵ0)---(20)]]>as=σsks(σs+ksσs)(bs-1)---(21)]]>σs=σsmax|s-s0|mdm---(22)]]>ks=1+(ksmax-1)|s-s0|mdm---(23)]]>αs=αsmax|s-s0|mαdmα---(24)]]>σmax=σopt≈m+1150πΔ---(25)]]>kzmax＝8  (26)kxmax＝kymax＝1  (27)m＝4  (28)其中m与mα是多项式的阶数，s0是CPML与模拟区域交界位置，s是当前CPML位置，d是CPML层的厚度，σ、α为中间计算变量，用于as、bs的计算；根据以上计算式分别得到各个场分量计算方程的系数和方程右端常数项，根据计算方程的特点，其系数构成了三对角矩阵，将三对角矩阵压缩存储在GPU全局内存中的数组A、B、C中，方程右端常数项存储在GPU全局内存中的数组X中；步骤3、使用在GPU上开辟的线程去求解计算电场的矩阵方程，从而得到每个网格处的电场场值Ex、Ey、Ez，最后进行电场边界的处理；步骤4、利用步骤3计算得到的电场场值，使用在GPU上开辟的线程去求解计算磁场的矩阵方程，从而得到每个网格处的磁场场值Hx、Hy、Hz，最后进行磁场边界的处理；步骤5、存储计算结果至GPU全局内存并重复步骤3、4，在每个时间步长都更新每个网格的电场与磁场的场值，从step＝0迭代直至达到step＝STEPS，STEPS为设置的模拟总时间步长数；步骤6、将步骤5的计算结果从GPU传回CPU主机端，进行数据分析与存储。