非结构网格瞬态电磁场计算中的高效通信方法

李岷轩，江树刚，吴庆恺，林中朝

(西安电子科技大学 陕西省超大规模电磁计算重点实验室，陕西 西安 710071)

随着无线电技术和微电子技术的飞速发展，大型设备和电子电路的电磁兼容问题变得越来越普遍且重要。以电磁场仿真算法为基础的电磁场仿真软件作为天线设计、电磁兼容等工程问题的基础工具，具有非常重要的实际意义。在计算电磁学领域的众多电磁算法中，时域有限差分(Finite Difference Time Domain，FDTD)方法、有限元(Finite Element Method，FEM)方法和矩量法(Method of Moments，MoM)是3种精度较高的经典算法。而时域间断伽略金(Discontinuous Galerkin Time Domain，DGTD)方法兼具时域有限差分方法的时域宽频带特性和时域有限元(Finite Element Time Domain，FETD)方法的非结构网格特性，具有很高的理论精度，尤其在计算宽带场路耦合问题和复杂媒质问题时更具优势[1-7]。

不同于时域有限元方法求解强耦合单元的特性，时域间断伽略金方法使用间断伽略金法降低了网格间的耦合性，求解过程中计算各单元的独立矩阵，避免了计算大型稀疏矩阵，计算复杂度远小于时域有限元方法[8]。时域间断伽略金方法中常用的显式迭代方案有leapfrog和ADI等[9-10]，在时域间断伽略金方法的显式迭代并行计算过程中，各个单元在迭代更新时所需的相邻单元信息为前一时刻信息，因此需要与该单元的相邻单元进行信息交换。在并行迭代计算过程中，这一信息交换过程是时域间断伽略金方法的主要通信热点。相邻单元的数量会明显影响并行计算过程中的通信量，即通信负载。将单元划分给各并行进程的过程决定了各个进程的计算量，即各进程的计算负载。通信负载与计算负载共同决定了算法的并行效率。受制于区域划分方法的局限性[11-12]，与采用结构性网格进行离散的时域有限差分[13-14]方法相比，在时域间断伽略金方法并行计算的过程中需要在计算负载和通信负载间做取舍[15-16]。同时在并行迭代过程中，进程间的数据收发过程非常频繁，因此在大规模并行状态下难以取得很高的并行效率[15]。笔者将给出一种适用于大规模并行的时域间断伽略金通信策略。

1 并行DGTD原理与通信策略

无源麦克斯韦方程组的微分形式中，电场E和磁场H表示为

(1)

(2)

(3)

其中，Z为单元i的阻抗，Y为单元i的导纳；Z+为单元i相邻单元的阻抗，Y+为单元i相邻单元的导纳；E+、H+为单元i相邻单元的电场和磁场。由于数值通量的存在，求解式(2)的过程中需要相邻单元的电磁场值。对式(2)进行伽略金检验后得到[9]：

(4)

(5)

其中，Γlp为与单元i属于相同进程的相邻单元面，Γcp为与单元i属于不同进程的相邻单元面。以式(5)为基础，采用不同时间差分近似方法可以得到关于ei和hi的并行迭代方程。同时，在不同进程的交界面单元处，各进程的边界单元需要与相邻进程单元交换信息。

(6)

除边缘进程外(被分配区域处于模型边缘的进程)，各进程的通信量Pmessage与计算量Pamount之比为

(7)

其中，αcpn为相邻进程数，对于立方体均匀分割的情况，αcpn可视为6。γ为每个相邻网格所需的通信量，β为各个网格每次迭代时的计算量。式(7)表明，对每个进程来说，其通信量与进程数正相关，计算量与进程数负相关，当计算规模(进程数)扩大时，程序的计算效率会因为通信量的提升而降低。算法每次迭代所消耗的时间满足

(8)

其中，T为每次迭代的总时间，Tm为每次迭代的通信时间，Tc为每次迭代的计算时间，Sm为计算平台每个通信周期的通信速度，γp为各进程在每个通信周期中的平均通信量，Sc为计算平台的运算速度，N为整个通信过程的非自锁点对点通信次数(通信周期)。显然，要减少算法的迭代时间，需要提高算法的计算速度和通信速度。使用图形处理器(Graphic Processor Unit，GPU)和多核处理器(Many Integerated Core，MIC)卡等异构加速技术和中央处理器(Central Processing Unit，CPU)多线程技术能够直接而有效地提升算法的计算速度[15-17]；笔者着重讨论提升算法通信速度的方式。根据文献[10]中并行效率的定义，n个进程计算时的并行效率Rn可以表示为

(9)

其中，Tc为基准核数的计算时间，c为基准核数。视Sm、γp、Sc、β为常数，可以得到：

(10)

(11)

(12)

式(10)表明，在均匀划分的理想情况下，点对点通信的DGTD方法迭代过程的并行效率与并行规模和通信周期呈负相关。其中，Cm是一次迭代中方形结构一个面上所有网格的通信时间，Cc是单进程情况下整个模型的总迭代计算时间。显然，在大规模并行时想要取得较好的并行效率，就需要减小每次通信的通信周期。

基于消息传递接口(Message Passing Interface，MPI)中的MPI_ISEND和MPI_IRECV两个接口，可以方便地完成相邻进程的信息交换[12]。为了降低算法在通信时的通信周期，需要在通信过程开始之前规定好各进程点对点通信的顺序。最简单的通信顺序为进程号由小到大依次通信的顺序通信方式。以图1(a)所示的进程关联情况为例，其进程间的相邻关系可以用图1(b)所示的关联矩阵表示。图1(b)中列号为本地进程号，行号为相邻进程号，进程间有相邻关系的元素为1，无相邻关系的元素为0。采用顺序通信方式后，各个进程的通信顺序即通信拓扑如图1(c)所示，各个矩阵元素值代表对应进程之间所属的通信周期数；该策略通信过程中的实际通信情况如图2所示，图中由虚线箭头连接的方框表示一组收发过程，横轴为通信周期数。

实际应用过程中，一个进程往往与多个进程相关联，顺序通信方式下总通信周期数近似为进程总数，每个通信周期内的通信量远小于通信缓冲区大小，造成了通信资源的浪费。

当一对相邻进程进行通信时，与这两个进程没有相邻关系的另一对相邻进程称为这两个进程的无关通信组。在同一个通信周期内，几个无关通信组能够同时进行通信操作，从而避免通信带宽的浪费。如图3(a)所示，按进程号从小到大的顺序，依次分析每个通信周期内的无关通信组，能够得到通信周期最小的通信拓扑，称为最小周期通信拓扑。以下给出获取最小周期通信拓扑的伪代码。

最小周期通信拓扑伪代码：

DO WHILE NotAllCommunicationFinished

NowPeriod++；

DO send_rank=0，mpi_size

IF comm_check(send_rank).equal..false. THEN

DO recv_rank=ConnectProcessor0，ConnectProcessorN；

IF comm_check(recv_rank).equal..false. THEN

comm_check(send_rank)=.true.；CommMatrix(send_rank，recv_rank)=NowPeriod；

comm_check(recv_rank)=.true.；CommMatrix(recv_rank，send_rank)=NowPeriod；exit；

ENDIF

ENDDO

ENDIF

ENDDO

ENDDO。

最小周期通信的实际通信过程如图3(b)所示。相对于点对点通信中简单的顺序通信方式，最小周期通信方式的总通信周期更少；根据式(10)，其理论并行效率也更高。

2 数值算例

以图4(a)所示的矩形波导缝隙天线辐射特性计算为例，验证方法的有效性。矩形波导尺寸为710 mm×13.2 mm×6.5 mm，使用同轴线在波导中心处馈电，采用频带为14～20 GHz的调制高斯脉冲激励源激励，计算时使用由Intel Xeon Gold 6140 CPU @ 2.30 GHz处理器组成的计算集群；该集群共16节点，每节点含两块CPU处理器及72 GB内存；每个处理器核心数为18，最大线程数为36。计算时，给每个处理器核心分配两个进程。其中，缝隙天线16 GHz的xOz面和yOz面增益如图4(b)和图4(c)所示。表1给出了使用顺序通信和最小周期通信的计算情况。可以看出，在不同计算核心数并行规模的情况下，顺序通信的通信周期远大于最小周期通信，导致并行效率大大降低，这一点在大规模计算时更为突出。

表1 矩形波导缝隙天线计算情况(迭代663 217步)

文献[15]采用非阻塞的顺序通信策略，在计算1.6×106网格量的模型时，由6进程扩展至30进程(5倍)时效率降至90%以下。文中算例网格量为1.1×106，与文献[15]中算例有一定相似性。文中算例在核心数为936的计算过程中，每个进程平均计算的网格量为1 176，每对进程的单次通信的网格量约为121个单元的信息。以120进程为基准，采用最小周期通信策略扩展8倍的并行效率达到了约89%，同时其通信周期约为顺序通信的3%，这表明，采用最小周期通信策略能有效降低程序的通信周期，从而提升算法的并行效率。

以图5(a)所示的汽车模型为例，验证并行规模较大时，最小周期通信策略的并行效率。汽车尺寸为4.8 m×1.8 m×1.45 m，计算频率为1 GHz，计算平台为神威·太湖之光超级计算机，激励源为水平极化平面波，入射方向沿z轴负向，极化方向沿x轴正向，计算θ=90°，φ由0°到360°的xOy面双站散射截面(Radar Cross Section，RCS)。需要注意的是，神威·太湖之光超级计算机采用由申威26010异构众核处理器组成的计算节点，该处理器由1个主核和64个从核组成，称为一个核组，每个节点包含4个核组[18]。计算过程中，每个核组作为一个进程共同参与计算。DGTD方法在申威26010异构众核处理器上运行时需要对主核和从核分别进行编程[13]，而笔者主要讨论核组间的通信过程，对此不多叙述。

图5(b)展示了DGTD与MoM的RCS计算结果。该算例迭代10 000步的计算时间和并行效率如图6所示。以800进程为基准，当并行规模扩展到8 000进程(并行规模扩展10倍)时，并行效率约为70.38%。根据式(8)，图6中的通信时间受通信周期数和进程数增加的共同影响而逐渐上升，这是影响其并行效率的主要因素。

3 结束语

在并行规模较大时，点对点通信的DGTD方法的并行效率主要取决于每次迭代时通信过程的通信周期数。笔者采用最小周期通信策略能够有效降低通信时的通信周期，使算法在进程间进行点对点通信时取得更好的并行效率，从而提升DGTD方法在大规模并行计算时的计算效率，进而使其更高效地应用于工程问题的计算中。