一种基于迭代方法的并行化时域混合电磁算法与流程

本发明涉及电磁场数值计算技术领域，尤其涉及一种基于迭代方法的并行化时域混合电磁算法。

背景技术：

雷达系统对识别目标和远程预测报警需要通过目标的雷达散射截面(rcs)来提供数据支撑和准确判断，大尺寸复杂形状目标rcs计算问题可通过时域电磁仿真进行一次计算获得宽频带仿真结果，并行化则有利于计算规模的扩大及计算时间的缩短。单纯使用并行化时域有限差分法(fdtd)计算电磁问题时，整个计算区域为一个长方体，网格均匀剖分为六面体网格，按虚拟拓扑方式将网格平均划分到各进程，以平衡各进程负载，使得网格交界面处的网格量达到最少即为最优并行方案。而单纯使用并行化有限体积法(fvtd)方法时，计算区域的网格剖分为四面体网格，按邻接区域平均分配网格到各个进程，并优化交界面网格数使其最少即可保证最优的并行加速比。

现有技术中，针对复杂电磁结构，单纯采用六面体网格剖分会将模型的边界阶梯化导致建模误差，单纯采用四面体网格剖分可实现精确建模但往往网格量巨大。采用混合网格建模计算则是一个有效手段，对于弯曲、细小或者倾斜结构采用四面体网格剖分，模型的其余部分采用六面体网格剖分，既能保证建模精度网格量又不至于过大。混合网格计算方法具有建模上的优势，但由于涉及到两种网格，两种计算方法，导致其并行化过程变得复杂。并行技术的使用会使得随着参与计算的进程数目不断增加计算时间迅速减少。而计算规模增大时，场值通信的交界面网格总量增加，从而取代通信机制效益成为影响并行效率的主要因素。而且，对于相同进程数目的情况下，并行方案不同时计算时间差别很大。因此，如何在有限的计算资源下获得最佳并行方案成为混合方法中需要解决的重要问题。

技术实现要素：

针对模型阶梯化剖分、网格量巨大、并行过程复杂和存在两种网格的两种方法的计算问题，本发明的目的在于提供一种解决这些问题的基于迭代方法的并行化时域混合电磁算法。

为了实现上述目标，本发明采用以下技术方案：

步骤1，根据目标的最高求解频率计算网格尺寸；

步骤2，使用封闭的正方形面片包围面包围目标表面；

步骤3，将计算区域分为两部分：采用四面体网格剖分包围面至目标表面的区域以实现精确建模，使用六面体网格剖分包围面至计算区域边界的部分以节省计算资源；

步骤4，根据网格上所采用的迭代方法及所消耗的计算量，以均衡负载为原则，将参与计算的节点分为不同的迭代类别，计算每个节点上的网格量；

步骤5，根据网格上正常迭代需要用到的相邻网格物理量估算节点间的数据通信量，以数据通信量最少为设计目标，将网格分配到具体节点上实现最优并行方案；

步骤6，采用并行fdtd-fvtd混合方法迭代求解，得出计算结果。

进一步地，所述的根据目标的最高求解频率计算网格尺寸，方法为：

对于目标在[fmin，fmax]频率范围内雷达散射截面rcs，采用以下公式计算fdtd方法的网格尺寸δ：

δ＝(c0/fmax)/10

其中，c0＝3.0×10⁸m/s为光速。

进一步地，所述的使用封闭的正方形面片包围面包围目标表面，包括：

根据网格尺寸的大小，对目标表面进行三角形网格剖分，循环三角形网格面片中心点坐标外推2-3个网格尺寸建立由边长为所述网格尺寸大小的正方形面片构成的包围面。

进一步地，将参与计算的节点分为不同的迭代类别，包括：

所有四面体网格仅进行fvtd迭代，与四面体网格相邻的六面体网格也仅进行fvtd迭代，与这些仅进行fvtd迭代的六面体网格有公共面的六面体网格称之为过渡区，这部分网格上既进行fdtd迭代又进行fvtd迭代，除此之外，过渡区外的六面体网格则仅进行fdtd迭代。

进一步地，所述的计算每个节点上的网格量，包括：

设仅进行fdtd方法迭代的六面体网格个数为n1，仅进行fvtd方法迭代的四面体网格个数为n2，仅进行fvtd方法迭代的六面体网格个数为n3，既进行fdtd迭代又进行fvtd迭代的六面体网格个数为n4，则：

fdtd迭代区域所有网格上的计算量为c1＝7×6×n1；

其余部分网格上的计算量c2＝97×6×n2+145×6×n3+(145+7)×6×n4。

进一步地，步骤4还包括：将所有网格按照估算的计算量分配节点：

为fdtd迭代区域分配d个节点，其计算公式为：

d＝[c1×m/(c1+c2)]

为fvtd迭代区域分配v个节点，其计算公式为：

v＝[c2×m/(c1+c2)]

d和v之间满足关系：d+v＝m，m为节点总数；

将n1个参与fdtd迭代的六面体网格平均分布到第0-d-1号节点上；将(n2+n3+n4)个参加fvtd迭代的网格平均分布到第d—m-1号节点上。

进一步地，所述的根据网格上正常迭代需要用到的相邻网格物理量估算节点间的数据通信量，包括：

节点i(0≤i≤d-1)和其相邻节点j(0≤j≤d-1,且i≠j)上分布的网格均进行fdtd方法迭代，设两个节点公共边界上有相邻网格n1个，该部分网格的数据通信量为：2×n1×4；

节点i(0≤i≤d-1)上分布的网格进行fdtd方法迭代，若其相邻节点j满足(d≤j≤m-1,且i≠j),即节点j上的网格进行fvtd方法迭代，设两个节点公共边界上有相邻网格有n2个，该部分网格的数据通信量为：2×n2×4；

节点i(d≤i≤m-1)上分布的网格进行fvtd方法迭代，相邻节点j(d≤j≤m-1,且i≠j)上的网格也进行fvtd方法迭代，设两个节点公共边界上有相邻网格有n3个，该部分网格的数据通信量为：2×n3×26。

进一步地，所述的以数据通信量最少为设计目标，将网格分配到具体节点上实现最优并行方案，包括：

将0—d-1号节点进行三维虚拟拓扑建模，即设x方向分布a个节点，y方向分布b个节点，z方向分布c个节点，a,b,c均为整数，且由以下两个条件确定：

a:b:c≈(lxmax-lxmin):(lymax-lymin):(lzmax-lzmin)

a×b×c＝d

节点的三维拓扑坐标为：(0,0,0)，(1,0,0)…(a-1,b-1,c-1)；

将中心点坐标(x,y,z)满足条件：

lxmin+(o-1)δ≤x≤lxmin+oδ

lymin+(p-1)δ≤y≤lymin+pδ

lzmin+(q-1)δ≤z≤lzmin+qδ的网格预分布到拓扑坐标为(o,p,q)的节点上；

从第(0,0,0)拓扑节点开始，遍历每一个节点执行以下步骤：

若节点上预分布的网格个数大于n1/d，则将该节点边界面上的网格全部分配到相邻节点上；

若节点上预分布的网格个数小于n1/d，则将相邻节点边界面上的网格全部分配到该节点上；

按照网格编号顺序由编号小的网格按照邻接关系向外围遍历的方式分布到节点d—m-1上；

估算不同拓扑结构下节点间的数据通信量，选择数据通信量最小的分配方案为最终的并行方案。

本发明与现有技术相比具有以下技术效果：

1.本发明提出一种并行化的电磁场计算方法，避免了阶梯化剖分以及网格量大的问题，并行化方案最大化利用了计算机资源缩短了计算时间。

2.本发明根据具体模型的网格分布特点及网格上所进行的迭代方法不同对计算区域进行划分，以进程负载平衡为设计准则，以减少进程间的数据通信量提高并行加速比为设计目标，为大规模复杂几何模型电磁问题的求解提供一种简单高效的并行化方案。

附图说明

图1为球形目标的正方形面片包围面示意图；

图2为采用四面体网格剖分包围面至球形目标表面之间的区域剖面图；

图3为采用六面体网格剖分包围面至计算区域边界之间的区域剖面图；

图4为整个计算区域的网格剖分及网格上所采用的迭代方法示意图；

图5为包围面附近网格的迭代方法示意图；

图6为按网格上所采用迭代方法划分计算区域的示意图；

图7为按网格邻接关系划分计算区域的示意图；

图8为球形目标的双站rcs计算结果；

图9为本发明方法的流程示意图。

具体实施方式

本发明对存在弯曲，精细或倾斜结构的大规模电磁场数值仿真问题，使用并行化时域有限差分法和有限体积法(fdtd-fvtd)的混合方法求解时涉及两种网格两种迭代方式数值方法的并行方案设计。应用广泛的并行fdtd方法将计算区域剖分为六面体网格，按照虚拟拓扑方式设计的并行方案简单直接效果好，但对于存在两种网格的两种方法的计算问题，并行方案的设计研究不足。本发明基于网格所进行的迭代方法不同对计算区域进行划分，以进程负载平衡为设计目标，尽量减少进程间的数据通信量，提高并行加速比，进而达到在计算前使用该方法可获得在有限计算资源下适用于fdtd-fvtd混合方法计算的最优并行方案。

如图1至图8所示，本发明公开了一种基于迭代方法的并行化时域混合电磁算法，以计算一个目标在[fmin，fmax]频率范围内rcs结果为例，由于计算机容量的限制，数值计算只能在有限区域进行。为了模拟开域电磁问题，在计算区域的截断边界处必须采用吸收边界。吸收边界条件所确定的问题空间是有限的，经合理设置(本发明中的混合方法沿用了fdtd方法中经典的完美匹配层pml)后的有限空间可与无限空间等效；这个有限空间即为计算区域，它是沿坐标轴x,y,z三个方向建立，包含目标在内，边界由pml构成的一个长方体空间，设其在三个方向的顶点上下限分别为[lxmin,lxmax]，[lymin,lymax]，[lzmin,lzmax]。该方法的详细步骤如下：

步骤1，根据目标的最高求解频率计算网格尺寸；

对于目标在[fmin，fmax]频率范围内雷达散射截面rcs，采用以下公式计算fdtd方法的网格尺寸δ：

δ＝(c0/fmax)/10

其中，c0＝3.0×10⁸m/s为光速。

步骤2，使用封闭的正方形面片包围面包围目标表面；

混合方法的实现需要将整个计算区域剖分为两种网格，计算区域中哪一部分使用四面体网格剖分，哪一部分使用六面体网格剖分需要在剖分前确定好。因此，这一步骤的目的是将计算区域分成两部分，用于后面两种网格的剖分。

按照步骤1获得的网格尺寸大小(边长)δ，使用网格剖分软件对目标表面进行三角形剖分，循环三角形网格面片中心点坐标(以此方式近似获得目标表面的几何坐标)外推2δ-3δ距离，建立由边长为δ的正方形构成的包围面。该包围面将作为下一步骤中网格剖分时四面体网格和六面体网格的交界面。图1所示为球形目标模型的包围面。

步骤3，将计算区域分为两部分：采用四面体网格剖分包围面至目标表面的区域以实现精确建模，使用六面体网格剖分包围面至计算区域边界的部分以节省计算资源；

利用上一步骤中获得的包围面将整个计算区域分成两个部分：包围面至目标表面和包围面至计算区域边界。

步骤3.1，使用网格剖分软件将包围面与目标表面之间的区域进行四面体网格剖分，四面体网格可精确描述目标的几何形状，包围面和目标间的距离仅2-3个δ限制了四面体网格的大范围使用，进而使得网格量不至于过多，这是混合方法的一个优点。如图2所示为包围面至球形目标表面之间的区域被剖分为四面体网格的剖面图。

步骤3.2，使用网格剖分软件将包围面外围至计算区域边界的部分采用六面体网格剖分。如图3所示为球形目标包围面外围至计算区域边界区域被剖分为六面体网格的剖面图。混合方法的另一个优点是计算区域内大范围使用的网格是六面体网格采用简单直接的fdtd方法迭代，兼顾了计算效率。

至此，整个计算区域的网格剖分完成，对六面体网格进行一维编号(1,2，…)的同时计算其三维编号，如某网格中心的坐标为(x,y,z)，则该网格的三维编号为(i,j,k)，具体计算公式为：

其中代表对其内部变量进行向上取整。如图4所示为球形目标的网格剖分的剖面示意图。

步骤4，根据网格上所采用的迭代方法及所消耗的计算量，以均衡负载为原则，将参与计算的节点分成fdtd方法迭代和fvtd方法迭代两类，计算每个节点上的网格量；

步骤4.1，fdtd和fvtd两种方法均揭示了场的局部特性，二者的混合只需要在两种方法网格交界处进行场值交换即可实现。如图5所示，所有四面体网格仅进行fvtd迭代，与四面体网格相邻的六面体网格也仅进行fvtd迭代，与这些仅进行fvtd迭代的六面体网格有公共面的六面体网格称之为过渡区，这部分网格上既进行fdtd迭代又进行fvtd迭代，除此之外，过渡区外的六面体网格(包括pml层)则仅进行fdtd迭代。球形目标rcs计算问题中，整个计算区域中各个网格上所采用的迭代方法如图4所示。

步骤4.2，节点上网格量的计算

设仅进行fdtd方法迭代的六面体网格个数为n1，仅进行fvtd方法迭代的四面体网格个数为n2，仅进行fvtd方法迭代的六面体网格个数为n3，既进行fdtd迭代又进行fvtd迭代的六面体网格个数为n4。

fdtd方法中，一个六面体网格单元上需要6个场值(电场的三个分量：ex,ey,ez以及磁场的三个分量：hx,hy,hz)，实现一次迭代需要进行的计算次数大约为7×6。fvtd方法中，采用二阶迎风格式进行重构时，考虑单元间的邻接关系，一个四面体单元实现一次迭代需要的平均计算次数为：97×6，一个六面体单元实现一次迭代需要的平均计算次数为145×6。

因此，设fdtd迭代区域所有网格上的计算量为c1，其计算公式为：

c1＝7×6×n1

其余部分网格计算量c2的计算公式为：

c2＝97×6×n2+145×6×n3+(145+7)×6×n4

步骤4.3，设参与并行计算的节点个数为m(编号分别为：0,1,2…m-1)，考虑网格间的邻接关系及网格上所采用的迭代方法，将所有网格按照估算的计算量分配节点：

设为fdtd迭代区域分配d个节点，其计算公式为：

d＝[c1×m/(c1+c2)]

设为fvtd迭代区域分配v个节点，其计算公式为：

v＝[c2×m/(c1+c2)]

d和v之间满足关系：d+v＝m，其中[]代表对其内部变量进行四舍五入取整。

将n1个参与fdtd迭代的六面体网格平均分布到第0-d-1号节点上，每个节点上分布约n1/d个网格；将(n2+n3+n4)个参加fvtd迭代的网格平均分布到第d—m-1号节点上，每个节点上分布约(n2+n3+n4)/v个网格。

该步骤确定了每个节点上分布的网格个数，可使计算量均匀分布到各个节点上，达到负载均衡的效果。

步骤5，根据网格上正常迭代需要用到的相邻网格物理量估算节点间的数据通信量，以数据通信量最少为设计目标，将网格分配到具体节点上实现最优并行方案。

存在公共面的体网格称之为相邻网格。一个网格正常迭代时需要用到相邻网格上的场量及其他物理量，节点交界面上的网格迭代时需要其相邻网格(分布在该节点的相邻节点上)的场值和物理量，为使迭代正常进行，需要将这些场值和物理量通过通信方式发送。数据通信量越多，通信时间占计算时间的比例就越大，严重影响并行加速比，因此，在计算前估算数据通信量将利于最优并行方案的设计。

fdtd方法中位于节点交界面上的每个网格需要发送2个场值数据，同时接收2个场值数据。fvtd方法中位于节点交界面处的每个网格需要发送13个值，同时接收13个值。采用fvtd方法迭代的网格交界面场值通信量远大于fdtd方法迭代的网格交界面场值通信量。因此，减少fvtd方法迭代区域的子节点划分数量可大大降低通信时间的消耗。

步骤5.1，数据通信量的计算分以下三种情况：

①节点i(0≤i≤d-1)和其相邻节点j(0≤j≤d-1,且i≠j)上分布的网格均进行fdtd方法迭代，设两个节点公共边界上有相邻网格n1个，该部分网格的数据通信量为：2×n1×4。

②节点i(0≤i≤d-1)上分布的网格进行fdtd方法迭代，若其相邻节点j满足(d≤j≤m-1,且i≠j)，即节点j上的网格进行fvtd方法迭代，设两个节点公共边界上有相邻网格有n2个，该部分网格的数据通信量为：2×n2×4。

③节点i(d≤i≤m-1)上分布的网格进行fvtd方法迭代，相邻节点j(d≤j≤m-1,且i≠j)上的网格也进行fvtd方法迭代，设两个节点公共边界上有相邻网格有n3个，该部分网格的数据通信量为：2×n3×26。

下一步骤中并行方案的设计将以节点交界面上所有网格的数据通信量总和最少作为设计目标。

步骤5.2，并行方案的设计

根据步骤4.2计算得到的各节点上网格个数，以场值通信量总和最少为目标将每个网格分配到具体节点上是并行方案设计的最后一步。具体分配方案如下：

步骤5.2.1，fdtd迭代区域网格的分配方案：

将0—d-1号节点进行三维虚拟拓扑建模，即设x方向分布a个节点，y方向分布b个节点，z方向分布c个节点，a,b,c均为整数，且由以下两个条件确定：

a:b:c≈(lxmax-lxmin):(lymax-lymin):(lzmax-lzmin)

a×b×c＝d

节点的三维拓扑坐标为：(0,0,0)，(1,0,0)…(a-1,b-1,c-1)。

将中心点坐标(x,y,z)满足条件：

lxmin+(o-1)δ≤x≤lxmin+oδ

lymin+(p-1)δ≤y≤lymin+pδ

lzmin+(q-1)δ≤z≤lzmin+qδ的网格预分布到拓扑坐标为(o,p,q)的节点上；

如图6所示，采用6个计算节点计算球形目标双站rcs时按照x方向分布3个计算节点，y方向分布2个计算节点，z方向分布1个计算节点的网格分布及各个节点的拓扑坐标。

考虑节点上分布的网格个数要均衡，由步骤4.3可知，每个节点上将分布n1/d个网格。由于计算区域内并不完全是整齐排列的六面体网格，还有一部分四面体网格，因此按拓扑结构分布的各节点上网格个数并不均匀，还需要调整，以使各节点上的计算量均衡：

步骤5.2.2，从第(0,0,0)拓扑节点开始，遍历每一个节点执行以下步骤：

若节点上预分布的网格个数大于n1/d，则将该节点边界面上的网格全部分配到相邻节点上；

若节点上预分布的网格个数小于n1/d，则将相邻节点边界面上的网格全部分配到该节点上。

步骤5.2.3，由步4.3可知，节点i(d≤i≤m-1)上将分布(n2+n3+n4)/v个fvtd迭代区域的网格，按照网格编号顺序由编号小的网格按照邻接关系向外围遍历的方式分布到节点d—m-1上；

最后，按照步骤5.1的数据通信量计算方式估算不同拓扑结构(即不同a,b,c组合)下节点间的数据通信总量，选择数据通信量最小的分配方案为最终的并行方案。

针对球形目标的rcs计算问题，如采用6个节点的并行fdtd-fvtd混合方法计算时，如图7所示为使用本发明提出的并行方案得到的计算区域划分剖面图。根据网格所采用的迭代方法将四面体网格区域和六面体网格区域分别划分，本方法考虑了四面体网格和六面体网格对计算资源消耗的不同以及数据通信量对并行效率的影响。

如图6所示，fvtd迭代区域内部，不同节点的网格交界面有3个。图6为按照网格分布，单纯根据邻接关系不考虑网格所采用的迭代方法，将所有网格平均分布到所有节点上的网格划分情况，fvtd迭代区域内部不同节点的网格交界面有6个。由步骤5.2第3种情况可知，相比于前两种情况，fvtd迭代区域节点交界面的增加会大大增加数据通信量，影响并行效率。

步骤6，采用并行fdtd-fvtd混合方法迭代求解，得出计算结果。

使用步骤5得到的并行方案将网格分布到m个计算节点上，然后采用并行fdtd-fvtd混合方法迭代求解。所述的并行fdtd-fvtd混合方法的出处为：贺之莉，基于gid的并行fdtd及其混合方法研究，博士论文，2012年6月。

如图8为半径为1m的球形目标在频率为1ghz的双站rcs计算结果与矩量法(mom)计算结果吻合良好。