CN107526887B - 一种基于GPU并行的LeapfrogADI-FDTD方法 - Google Patents

Abstract

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

Description

一种基于GPU并行的LeapfrogADI-FDTD方法

技术领域

本发明属于计算电磁学与高性能计算领域。涉及一种基于GPU并行加速的计算技术，具体为一种基于GPU并行的Leapfrog ADI-FDTD方法。

背景技术

近30年来，由于高功率微波和毫米波源技术在科学技术和工业方面的广泛应用，高功率微波和毫米波源技术发展很快。近20年来，由于在生命医学、天体物理、等离子体加热、雷达和高速数据通讯等方面具有重要的应用前景，太赫兹技术也得到了快速的发展。无论是高功率微波和毫米波的应用，还是太赫兹的应用，中等功率或高功率的源是基础。真空电子器件是能够工作在室温条件下，并产生高功率微波、毫米波和低频段太赫兹波的主要器件。真空电子器件是通过将高压脉冲所具有的能量转化为电磁波的能量，通过波导以及天线将电磁波能量辐射出去。真空电子器件主要包括：返波管、回旋管、速调管、行波管、磁控管、振荡器和自由电子激光器等。

因为真空电子器件的实验设计和研制费用很高，研制周期较长，所以，现阶段真空电子器件的设计越来越依赖数值模拟，这样可以大大减少器件设计所需的时间及经费。在真空电子器件设计中的关键物理问题之一是微波的传输与放大特性，如何精确求解时变电磁场至关重要。在微波传输与放大特性的数值模拟中主要解决三个问题：1、内部场的求解问题；2、吸收边界的处理问题；3、计算效率问题。在真空电子器件内部，时变电磁场的求解采用的主要算法有时域有限差分算法(FDTD)、交替隐式时域有限差分算法(ADI-FDTD)以及蛙跳交替隐式时域有限差分算法(Leapfrog ADI-FDTD)。

FDTD是由Maxwell方程组离散后得到的显式差分算法，必须满足CFL时间稳定性条件，也就是时间步长的取值由计算空间中的空间步长的最小值决定，这就要求时间步长必须取的足够小。而时间步长的减小，使得迭代步数和运算时间大量增加，极大地增加了计算成本。为了弥补FDTD的缺陷，90年代后期曾提出了交替隐式时域有限差分算法(ADI-FDTD)。该算法存在中间的子时间步，在无条件稳定的同时损失了一定的计算资源，同时也使得对其并行变得无法实现，制约了它的发展与应用范围。

针对ADI-FDTD算法出现的不足，在2009年Cooke等人提出了一种新的无条件稳定的蛙跳交替方向隐式时域有限差分(Leapfrog ADI-FDTD)算法。在Leapfrog ADI-FDTD算法中，E和H在时间格点上同样相差Δt/2，但E和H都采用隐式公式。其每个时间步长所计算的场量与FDTD算法和ADI-FDTD算法的对比图如图1所示。较传统的ADI-FDTD算法，该算法没有子时间步的计算资源消耗，迭代方程更为简单，但仍然保持了无条件稳定的优点，所以该算法相对传统的无条件稳定算法有更高的效率。

由于计算机的局限性，基于有限差分的电磁数值模拟只能在有限区域内进行，这样必须对计算区域的边界进行截断处理。分析现有文献发现，在数值模拟电磁场的过程中，通常假设计算区域足够大，使得截断边界离激励源和目标的距离足够远。该方式由于计算域范围过大，数值计算时占用大量的计算内存，耗费大量的计算时间，导致计算效率极低。

为了克服以上缺点，在使用FDTD算法模拟电磁场时，通常采用缩小计算域范围，并在截断边界处设置合适吸收边界的形式。目前应用较多的吸收边界有Mur吸收边界、廖氏吸收边界和完全匹配层(Perfectly Matched Layer，PML)吸收边界。

Mur吸收边界虽然简单方便、占用内存小，但计算精度不高，在角区域存在较大误差；二阶近似的Mur吸收边界虽然精度有所改善，但对低频感应场的吸收效果差。廖氏吸收边界能够很好地吸收来自计算区域的电磁波，但它的高阶吸收边界不稳定。PML吸收边界给吸收边界的研究带来了本质性的突破，但它采用分裂电磁场的方式来描述各向异性材料，对计算机的硬件要求很高，并且在模拟电磁波时对高有耗介质中的波和倏逝波的吸收效果不佳。而在PML基础上发展的UPML(uniaxial PML)吸收边界则直接使用各向异性材料而不改变电磁场形式，能作为电磁波的吸收边界，但它的递推公式与材料特性有关，并且对大角度入射波及电磁波里低频成份的吸收效果差。

目前最具潜力的当属卷积完全匹配层(Convolutional PML，CPML)吸收边界，它除了具有UPML吸收边界的性能外，还能对倏逝波、晚时反射波以及低频感应场进行有效的吸收，因而CPML吸收边界特别适合用做电磁波数值模拟计算中截断边界的处理。

简单来说，为了达到较好的数值模拟效果，目前最好的处理办法是将LeapfrogADI-FDTD算法与CPML吸收边界结合，形成统一的计算算法进行数值模拟计算，这样既保证了Leapfrog ADI-FDTD算法的优良计算特性，同时也兼顾了CPML的吸收特性。但是，由于这种隐式求解方式需要求解矩阵方程，而求解矩阵方程的过程中计算资源消耗大、计算效率低下、计算时间长，因此增加了数值模拟实验的时间周期和成本，无法达到相关科学研究或器件设计的高效需求。

发明内容

针对上述存在的问题和不足，为解决含有CPML吸收边界的Leapfrog ADI-FDTD算法计算效率低下、计算时间过长的问题，本发明提供了一种基于GPU并行的Leapfrog ADI-FDTD方法，其基于GPU并行且含有CPML吸收边界。具体技术方案如下：

步骤1、将求解区域使用Yee网格划分。

对CPU主机端与GPU上的数据进行初始化，即在CPU主机端定义结构体数组，大小为模拟求解区域的总网格数，用于存放GPU返回的各个网格点的电磁场计算结果，并将各个网格的电磁场场值赋值为0。同时，在GPU上的全局内存中定义同样大小的结构体数组，用于存储电磁场计算结果，并将各个网格的电磁场场值赋值为0。

步骤2、数据预处理，在GPU全局内存上开辟空间，并利用下述推导的计算式，计算得到矩阵方程的系数矩阵。数据预处理流程如图4所示。

在三维情形下，含有CPML的Leapfrog ADI-FDTD算法计算的场量有六个，包括：电场的三个分量：E_x、E_y、E_z；磁场的三个分量：H_x、H_y、H_z。其中，E_x的计算公式如式(1)、(2)、(3)所示，H_x的计算公式如式(4)、(5)、(6)所示：

其中：

其中:

其中，μ为磁导率，ε为介电常数，k_s为坐标伸缩系数(s＝x,y,z)，Δt为单位时间步长，Δx、Δy、Δz分别为沿x、y、z方向的单位网格长度，a_s与b_s为中间计算结果，二者的计算是为了得到ψ，其计算式如式(21)(20)所示，ψ为方程右端常数项的中间结果，i、j、k为模拟区域中网格在x、y、z方向上的编号，n为当前时间步长。其他方向的电场E_y、E_z与磁场分量H_y、H_z计算式同理可得。

由于在Yee网格中，存在半网格位置，为了在计算机中进行数值模拟计算，我们需要将地址映射到整数网格点。将上述(1)～(6)式进行地址映射与化简后我们得到E_x的计算式如式(7)～(12)所示，H_x的计算式如式(13)～(18)所示：

其中：

C_Ex＝A_Ex (10)

其中：

C_Hx＝A_Hx (16)

其中，A_Ex、B_Ex、C_Ex为E_x计算方程中左端的系数,A_Hx、B_Hx、C_Hx为H_x计算方程中左端的系数。其他方向电场分量E_y、E_z与磁场分量H_y、H_z计算式同理可得。

其中，方程中特定的变量求解式如下：

(s＝x,y,z) (19)

k_zmax＝8 (26)

k_xmax＝k_ymax＝1 (27)

m＝4 (28)

其中m与mα是多项式的阶数，s₀是CPML与模拟区域交界位置，s是当前CPML位置，d是CPML层的厚度，σ、α为中间计算变量，用于a_s、b_s的计算。

根据以上计算式分别得到各个场分量计算方程的系数和方程右端常数项，根据计算方程的特点，其系数构成了三对角矩阵，我们将三对角矩阵压缩存储在GPU全局内存中的数组A、B、C中，方程右端常数项存储在GPU全局内存中的数组X中。

步骤3、使用在GPU上开辟的线程去求解计算电场的矩阵方程。从而得到每个网格处的电场场值E_x、E_y、E_z，最后进行电场边界的处理。电场计算流程如图5所示。

步骤4、利用上一步计算得到的电场场值，使用在GPU上开辟的线程去求解计算磁场的矩阵方程，从而得到每个网格处的磁场场值H_x、H_y、H_z，最后进行磁场边界的处理。磁场计算流程如图6所示。

步骤5、存储计算结果至GPU全局内存并重复步骤3、4，在每个时间步长都更新每个网格的电场与磁场的场值，从step＝0迭代直至达到step＝STEPS，STEPS为设置的模拟总时间步长数。

步骤6、将步骤5的计算结果从GPU传回CPU主机端，进行数据分析与存储。

本发明适用于一维、二维及三维结构，适用于二维时，分布在激励源网格处的电磁场分量分别对应为E_y、H_z；电磁场分量对应公式相应退化。适用于一维时，分布在激励源网格处的电磁场分量对应为E_y；电磁场分量对应公式相应退化。

GPU是为了解决计算密集型、高度并行的问题而设计，它拥有大规模高度并行的流处理器，能够把高强度的计算需求通过合适的方式转变成为并行的流式数据计算模式，吞吐量很大，可以提供数十倍乃至于上百倍高于CPU的计算性能，在速度上获得一个数量级的提升。

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

附图说明

图1为三种算法中场量随时间步长更新方式对比图；

图2为三维Yee网格示意图；

图3为算法总体流程图；

图4为预处理模块流程图；

图5为电场求解模块流程图；

图6为磁场求解模块流程图；

图7为三维矩形波导计算实例示意图；

具体实施方式

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

实例中测试了三维均匀矩形波导模型。其横截面示意图如图7所示。波导宽度为a＝20.0mm，对模拟求解区域进行Yee网格划分，三个方向上的网格数分别为：N_x＝22、N_y＝12、N_z＝142，波导边界采用金属边界，取CPML层的厚度为40层(40*Δz)，加入CPML吸收边界的位置为沿传播方向Z的第100号网格，故令s₀＝99。激励源采用了TE₁₀模式的强迫性激励源作为面源，频率f＝10GHz，激励源加载在波导输入端口的截面上，所加激励源方程为：

ω＝2πf (30)

其中，X＝0.02，Δx＝1mm。单位网格大小Δx＝Δy＝Δz＝1mm，单位时间步长Δt＝0.57749ps，模拟时间步长数STEPS＝5000。

步骤1、将模拟求解区域使用Yee网格划分，三个方向上的网格数分别为：N_x＝22、N_y＝12、N_z＝142，在CPU主机端定义结构体数组，大小为22*12*142，用于存放GPU返回的电磁场计算结果并初始化各个网格的电磁场场值为0，同时在GPU的全局内存上定义同样大小的结构体数组，用于存储电磁场计算结果，并将各个网格的电磁场场值赋值为0。

步骤2、数据预处理，在GPU全局内存上开辟空间，根据前一部分推导的计算式我们可知，其系数构成了三对角矩阵，我们将三对角矩阵压缩存储在GPU全局内存中的数组A、B、C中，方程右端常数项存储在GPU全局内存中的数组X中。其中E_x的系数数组A、B、C与方程右端常数项中的ψ数组的计算式如(31)～(35)所示，H_x的系数数组A、B、C与方程右端常数项中的ψ数组的计算式如(36)～(40)所示：

C_Ex＝A_Ex (33)

C_Hx＝A_Hx (38)

其他方向的电磁场计算方程的系数数组A、B、C与常数项数组X也可同理得到。

步骤3、使用在GPU上开辟的线程去求解计算电场的矩阵方程。从而得到每个网格处的电场场值E_x、E_y、E_z，同时按金属边界条件将波导边界截面的切向电场分量赋值为0，最后根据当前时间步长更新激励源。

步骤4、利用步骤3计算得到的电场场值，使用在GPU上开辟的线程去求解计算磁场的矩阵方程，从而得到每个网格处的磁场场值H_x、H_y、H_z，最后按金属边界条件将波导边界截面的法向磁场分量赋值为0。

步骤5、将步骤3、4的计算结果存储在GPU的全局内存中并重复步骤3、4，从step＝0迭代直至达到STEPS＝5000。

步骤6、将步骤5的计算结果从GPU的全局内存传回CPU主机端内存中，保存计算结果。

经过计算并对比未进行GPU加速的CPU串行版本的含有CPML的Leapfrog ADI-FDTD算法可知，计算速度有了显著的提升。由CPU串行版本的16465.904秒计算时间缩短到GPU炳新版本的316.652秒计算时间，加速比达到了52，显著缩短了模拟计算时间，极大地提高了计算效率，减少了数值模拟实验的时间周期和成本，达到了相关科学研究或器件设计的高效需求。

Claims (3)

1.一种基于GPU并行的Leapfrog ADI-FDTD方法，具体步骤如下：

步骤1、将求解区域使用Yee网格划分；

对CPU主机端与GPU上的数据进行初始化，即在CPU主机端定义结构体数组，大小为模拟求解区域的总网格数，用于存放GPU返回的各个网格点的电磁场计算结果，并将各个网格的电磁场场值赋值为0；同时，在GPU上的全局内存中定义同样大小的结构体数组，用于存储电磁场计算结果，并将各个网格的电磁场场值赋值为0；

步骤2、数据预处理，在GPU全局内存上开辟空间，并利用下述推导的计算式，计算得到矩阵方程的系数矩阵；

含有CPML的Leapfrog ADI-FDTD算法计算的场量包括电场的三个分量E_x、E_y、E_z，磁场的三个分量H_x、H_y、H_z；其中，E_x的计算公式如式(1)、(2)、(3)所示，H_x的计算公式如式(4)、(5)、(6)所示：

其中，μ为磁导率，ε为介电常数，k_s为坐标伸缩系数s＝x,y,z，Δt为单位时间步长，Δx、Δy、Δz分别为沿x、y、z方向的单位网格长度，a_s与b_s为中间计算结果，二者的计算是为了得到ψ，其计算式如式(21)(20)所示，ψ为方程右端常数项的中间结果，i、j、k为模拟区域中网格在x、y、z方向上的编号，n为当前时间步长；

其他方向的电场E_y、E_z与磁场H_y、H_z分量计算式同理可得；

由于在Yee网格中，存在半网格位置，将地址映射到整数网格点，将(1)～(6)式进行地址映射与化简后得到E_x的计算式如式(7)～(12)所示，H_x的计算式如式(13)～(18)所示：

其中：

C_Ex＝A_Ex (10)

其中：

C_Hx＝A_Hx (16)

其中，A_Ex、B_Ex、C_Ex为E_x计算方程中左端的系数,A_Hx、B_Hx、C_Hx为H_x计算方程中左端的系数；其他方向电场分量与磁场分量计算式同理可得；

其中，方程中特定的变量求解式如下：

(s＝x,y,z) (19)

k_zmax＝8 (26)

k_xmax＝k_ymax＝1 (27)

m＝4 (28)

其中m与mα是多项式的阶数，s₀是CPML与模拟区域交界位置，s是当前CPML位置，d是CPML层的厚度，σ、α为中间计算变量，用于a_s、b_s的计算；

根据以上计算式分别得到各个场分量计算方程的系数和方程右端常数项，根据计算方程的特点，其系数构成了三对角矩阵，将三对角矩阵压缩存储在GPU全局内存中的数组A、B、C中，方程右端常数项存储在GPU全局内存中的数组X中；

步骤3、使用在GPU上开辟的线程去求解计算电场的矩阵方程，从而得到每个网格处的电场场值E_x、E_y、E_z，最后进行电场边界的处理；

步骤4、利用步骤3计算得到的电场场值，使用在GPU上开辟的线程去求解计算磁场的矩阵方程，从而得到每个网格处的磁场场值H_x、H_y、H_z，最后进行磁场边界的处理；

步骤5、存储计算结果至GPU全局内存并重复步骤3、4，在每个时间步长都更新每个网格的电场与磁场的场值，从step＝0迭代直至达到step＝STEPS，STEPS为设置的模拟总时间步长数；

步骤6、将步骤5的计算结果从GPU传回CPU主机端，进行数据分析与存储。

2.如权利要求1所述基于GPU并行的Leapfrog ADI-FDTD方法，其特征在于：适用于二维时，分布在激励源网格处的电磁场分量分别对应为E_y、H_z；电磁场分量对应公式相应退化。

3.如权利要求1所述基于GPU并行的Leapfrog ADI-FDTD方法，其特征在于：适用于一维时，分布在激励源网格处的电磁场分量对应为E_y；电磁场分量对应公式相应退化。

Images