一种针对复杂电磁环境的准三维抛物方程多级并行仿真方法

摘要

一种针对复杂电磁环境的准三维抛物方程多级并行仿真方法，采用三级并行方法实现对复杂电磁环境的快速仿真，即采用消息传递接口（MPI）进程级的任务划分实现外部二维切片的并行，再根据每个二维切片交叉步进的个数再进行相应的节点映射完成第二级并行，然后采用OpenMP线程级的任务划分与映射实现每个交叉步进上三对角方程组的并行，这样就把复杂的三维电磁环境仿真计算问题在简化成二维切片的粗粒度并行计算的基础上，再进一步对二维切片的计算过程完成了两级并行加速。相比简单采用粗粒度的MPI并行方法，三级并行方法使得计算速度更快，从而让并行计算的加速比和计算效率明显地提升，为接下来实现复杂电磁环境的更大范围和更快速度的仿真计算提供了有力基础。

说明书

技术领域

本发明涉及一种针对复杂电磁环境的准三维抛物方程多级并行仿真方法，即以电磁 传播与高性能计算等技术为支撑，在粗粒度并行计算的基础上更加快速的完成三维复杂 电磁环境仿真问题求解问题，例如微波与毫米波通信、雷达、精确制导、导航和地质勘 探等各种电磁领域，属于计算机领域。

背景技术

电磁环境是电磁干扰，电磁脉冲，电磁辐射对人员、军械和挥发性材料危害，以及 雷电和沉积静电等自然现象的综合。当处理大范围的复杂电磁环境仿真问题时，不管是 电磁干扰、电磁脉冲、电磁耦合的计算，复杂地形环境、大气、海洋等环境的建模，还 是现有的电磁仿真数值计算方法本身，都是一个极其复杂的过程，都需要耗费大量的时 间和资源，这些都对仿真过程中所用的计算机处理器、内存和存储器提出了很高要求， 而串行算法很难达到令人满意的效果，这就需要在复杂的计算中借助于高性能计算技 术，使仿真过程又快又好的完成，于是电磁环境高性能仿真应运而生，并且已经取得了 很多的成果。随着千万亿次超级计算机的研发成功，应用高性能计算技术来实现对复杂 电磁环境更快更好的仿真变得可能。电磁环境高性能仿真已经被广泛应用于诸如微波与 毫米波通信、雷达、精确制导、导航和地质勘探等各种电磁领域，具有巨大的实用价值。

目前基于抛物方程（Parabolic Equation，以下简称PE）的并行方法已有一些研究， 但是针对复杂电磁环境的并行方法，总体来说还是刚刚起步。PE方法在电磁环境仿真 领域越来越得到国内外学者的重视，是因为它和其它模型相比，本身就能反映电波传播 的折射和绕射效应。

PE是一种典型的偏微分方程，在数值求解过程中，有限差分法由于具有格式构造 简单、计算量少、应用灵活等突出优点，已成为最重要的数值方法之一。随着大规模科 学计算的需要和并行计算环境的发展成熟，数值分析工作者关于有限差分并行算法的研 究工作逐渐活跃起来。在交替分组显式(AGE)方法的基础上，张宝琳等提出利用Saul’yev 非对称格式构造分段隐式的思想，并恰当的使用交替技术建立了多种显-隐式和纯隐式交 替并行方法，取得了稳定性和并行兼顾的研究成果；在此基础上，吉林大学的吕桂霞等 针对二维变系数PE构造了一类变系数交替分块显-隐式（ABE-I）方法；山东大学李长 峰等结合迎风方法和区域分裂思想，采用一阶迎风、二阶修正迎风法逼近高维PE的对 流项、内边界处和子区域分别对应区域分裂显隐格式。在求解PE过程中出现了一些快 速迭代方法的研究，冯慧等通过不同点的隐式差分格式之间的相互约化来建立新型迭代 方法，此方法和Jacobi方法同样具有并行性，却比Jacobi收敛快；山东大学的张守慧等 在冯慧提出的数值Stencil的概念的基础上建立了收敛速度快、具有并行性的新型迭代格 式；近年来，还涌现出了大量的高阶交替分组格式的研究工作。

目前，并行PE方法现在已经应用在热传导以及其它扩散现象、化学反应、某些生物 形态、各种粒子的输运等领域，但是在电磁领域的应用才刚刚起步，需要国内外的研究 者对其展开广泛深入地研究。

针对准三维抛物方程求解大区域复杂电磁环境仿真过程中速度不够快的问题，寻找 一种复杂度较低的三维复杂电磁环境仿真的方法变得越来越迫切，本发明引入了多级并 行方法来加速求解的方法，该方法将每个二维切片的计算进行多级并行加速，采用交叉 步进并行进行第二级并行加速，再利用并行分裂法对三对角方程进行第三级并行加速。

发明内容

本发明技术解决问题：针对粗粒度的三维电磁环境并行仿真方法计算速度不够高、 计算区域不够大等问题，提供了一种针对复杂电磁环境的准三维抛物方程多级并行仿真 方法，在粗粒度并行计算的基础上更加快速的完成三维复杂电磁环境仿真问题求解问 题，这样就使并行计算效率进一步提高，最终在一定精度下满足实际生产需求中对于速 度的要求。

本发明提出一种针对复杂电磁环境的准三维抛物方程多级并行仿真方法，解决了微 波与毫米波通信、雷达、精确制导、导航和地质勘探等各种电磁领域的需求。本方法的 思想是在三维的二维切片划分和节点映射完成的基础上，根据每个二维切片交叉步进的 个数再进行相应的节点映射，在每个节点内部采用多核并行来加速求解三对角方程，把 计算的结果进行收集并进行相应的整理，最后显示出来。该方法主要包括四个步骤：二 维切片的划分及对二维切片与集群中的某几个节点进行相应映射的步骤；根据每个二维 切片交叉步进的个数再进行相应的节点映射的步骤；把每个节点内部采用多核并行来加 速求解三对角方程的步骤；计算的结果进行收集并进行相应整理的步骤。

具体包括以下：

步骤1，二维切片的划分及对二维切片与高性能计算集群中的某几个节点进行相应 映射，所述高性能计算集群是指64个计算节点，每个节点由2个Intel Xeon E5530处理 器和24GB的DDR31066GHz内存，总计算能力一共有512核处理器，完成对三维电磁 环境简化成二维切片的划分，并将划分好的二维切片与集群中的某几个节点进行相应的 映射，根据二维切片的个数和节点的特点进行对应的映射；

步骤2，在步骤1的二维切片划分和节点映射完成的基础上，根据每个二维切片交 叉步进的个数再进行相应的节点映射。根据计算过程所需的计算精度、集群的特点、二 维切片交叉步进的个数等因素，由天线的初始场来计算得出每个二维切片的初始值，将 相应的初始值分配给不同的计算节点；

步骤3，在步骤2交叉步进计算的节点映射完成的基础上，把每个节点内部采用多 核并行来加速求解三对角方程，在每个交叉步进的计算过程中，在每个节点上通过并行 分裂法对三对角线性方程组求解过程进行加速，在每个交叉步进计算对应的节点上对三 对角方程组的求解进行了加速；

步骤4，在步骤3的所有的线程并行计算完成的之后，计算的结果进行收集并进行 相应整理，根据每个二维切片在计算过程中划分成的很多交叉步进的计算，需要把计算 结果按照二维切片的交叉步进编号的进行交叉收集，然后按照每个二维切片排列成规则 的数据，最后实现数据可视化。

所述步骤2的根据每个二维切片交叉步进的个数再进行相应的节点映射，具体实现 过程如下：

在每个二维切片的计算过程中，首先根据计算过程中所需的计算精度来确定交叉步 进的个数，在这个过程中还需要考虑所用集群的节点、处理器和内存等特点，由初始场 的天线方向图来计算得出每个交叉步进上的初始值，采用如下公式进行计算： u(0,z)是初始场，θ是天线方向图的角度，F(θ)是天线方向 图函数，p是逆傅里叶变换的参数，H是天线高度，然后将相应的初始值分配给不同的 计算节点进行计算。

所述步骤3的在每个节点内部采用多核并行来加速求解三对角方程，具体实现过程 如下：

在每个交叉步进的计算过程中，在每个计算节点上需要求解一个三对角线性方程组 A_nU_n＝V_n＝A_n-1U_n-1，A_n是一个三对角矩阵，U_n是一个未知向量，V_n是一个已知向量，它由 前一个时间步的矩阵A_n-1和向量U_n-1相乘得到，通过并行分裂法对这个求解过程进行加 速，根据每个节点中的处理器个数对三对角矩阵进行等量的划分，然后通过在各个处理 器内的消去和处理器间的通信，将三对角矩阵变成对角形式，这样就能在各处理器内独 立的求解出方程的解，加速了三对角方程组的求解过程。

本发明与现有技术相比的优点在于：

（1）本发明把三维电磁环境空间进行划分后的二维切片进行了交叉步进并行的加 速，从而使得求解速度获得了接近线性的提高。

（2）本发明把每一个交叉步进上的三对角线性方程组的求解过程进一步采用并行 分裂法加速，从而让求解速度获得了进一步提升。

附图说明

图1为本发明基于抛物方程的准三维复杂环境多级并行方法框架；

图2为本发明中交叉步进求解方法示意图；

图3为本发明中采用分裂法在节点内的加速比示意图；

图4为本发明中抛物方程方法的加速比；

图5为本发明中抛物方程方法的并行效率；

图6为计算结果的显示效果图。

具体实施方式

下面结合附图对本发明作进一步详细的描述。

计算机集群（简称集群）是一种计算机系统，它通过一组松散集成的计算机软件和 /或硬件连接起来高度紧密地协作完成计算工作。在某种意义上，他们可以被看作是一台 计算机。集群系统中的单个计算机通常称为节点，通常通过局域网连接，但也有其它的 可能连接方式，单个节点内部可以有多个处理器组成。集群计算机通常用来改进单个计 算机的计算速度和/或可靠性。

天线方向图（也称为辐射方向图）是描述天线或其它信号源发出无线电波的强度与 方向（角度）之间依赖关系的图形。

本发明是一种针对复杂电磁环境的准三维抛物方程多级并行仿真算法，如图1所示， 该方法包括四个步骤：二维切片的划分及对二维切片与集群中的某几个节点进行相应映 射的步骤；根据每个二维切片交叉步进的个数再进行相应的节点映射的步骤；把每个节 点内部采用多核并行来加速求解三对角方程的步骤；计算的结果进行收集并进行相应整 理的步骤。

（1）二维切片的划分及对二维切片与集群中的某几个节点进行相应映射的步骤。

这个步骤是完成二维切片与集群节点的映射过程，如图1第一层并行的任务并行部 分，具体的实现过程如下所示：

①以发射机为中心，将基于抛物方程电磁环境传播三维空间的计算简化成一系列计 算相对简单的二维切片的划分；

②充分考虑所用的高性能计算集群的节点数、各节点的核数、计算能力、能耗等具 体特点，还需要根据计算精度的要求，完成计算过程中每个二维切片的交叉步进计算的 个数的划分；

③依据二维切片的编号和计算节点的编号，将划分好的各个二维切片与计算集群中 的某个或某几个计算节点进行相应的映射，从而完成第一层并行的任务并行部分。

（2）每个二维切片交叉步进的个数再进行相应的节点映射的步骤。

这个步骤是完成二维切片中的交叉步进计算与节点的映射过程，如图1第二层并行 的交叉步进并行部分，具体实现过程如下所示：

针对二维切片数远远大于计算节点数的情况，当各个二维切片在高性能计算集群上 的计算时，在计算过程所需的计算精度要求的范围之内，将这个计算过程分解成很多个 相对独立的二维切片交叉步进计算过程，这个相对独立的计算过程是将连续的时间步进 的过程分解成交叉的计算过程，每个二维切片的各交叉步进计算需要计算出相应的初始 场，计算过程如图2所示。在图2中可以看出，假如每个二维切片分解成5个交叉步进 进行分别计算。对于标准抛物方程：

$\frac{\partial u(x,z)}{\partial x}=\frac{\mathrm{ik}}{2}[\frac{2}{k^2}\frac{{\partial }^2}{\partial z^2}+n^2(x,z)-1]u(x,z)$

其中u(x,z)是在x0z平面的波函数，k是波数，n(x,z)是媒质折射系数。

根据计算的最小步长的大小，采用方向图函数的逆傅里叶变换，可以得到初始场为：

$u(0,z)=\sqrt{2\pi }{\int }_{-\infty }^{+\infty }F\left(\theta \right)e^{2\mathrm{\pi ipz}}\mathrm{dp}$

其中u(0,z)是初始场，θ是天线方向图的角度，F(θ)是天线方向图函数，p是逆傅里

变换的参数。

若将天线置于高度为H的地方，利用傅里叶变换中的平移定理，那么初始场为：

$u(0,z)=\sqrt{2\pi }{\int }_{-\infty }^{+\infty }F\left(\theta \right)e^{-2\mathrm{\pi ipH}}{e}^{2\mathrm{\pi ipz}}\mathrm{dp}$

其中H是天线高度，，计算出前4个初始场的值，将天线在内的5个初始场分配给 不同的计算节点，这样就使得每个二维切片的计算过程分解成了5个相对独立的计算过 程，计算速度大大提高，这个过程的理论基础是惠更斯原理。

（3）每个节点内部采用多核并行来加速求解三对角方程的步骤。

这个步骤是完成每个交叉步进计算中的方程并行计算过程，如图1第三层并行的方 程并行部分，具体实现过程如下所示：

在集群的每个计算节点的交叉步进并行计算过程中，需要完成对电磁环境的地形的 读取和插值、边界条件的计算、方程求解等工作，其中计算量最大的部分是复数的三对 角线性方程组的求解，一般耗费超过一半的计算时间，在矩阵的阶数越大的情况下就越 明显。相应的，集群中的各个计算节点上都有很多独立的核，可以进行采用OpenMP简 单的快速的完成对这个过程的加速计算，由此本发明采用比较精确的并行分裂法对这个 过程进行加速。并行分裂法就是把计算的矩阵分解成很多的块，然后将这些不同的块分 配到不同的处理器核上，通过核间的一些通信进行矩阵的消去操作，最后在每一个核上 都变成了一些对角方程的求解，这样就能方便的求解出三对角方程。

对标准抛物方程进行离散化和泰勒简化以后，可以得到如下三对角方程组的求解：

A_nU_n＝V_n＝A_n-1U_n-1

其中A_n是第n个时间步上的三对角矩阵；U_n是第n个时间步上的未知向量；V_n是第 n个时间步上的未知向量，它可由第n-1个时间步上的矩阵A_n-1和向量U_n-1的乘积得到， 具体形式如下所示：

$A_n\left(\begin{array}{cccccccc}1& 0& 0& 0& & ···& 0& 0\\ 1& {\alpha }_1^n& 1& 0& & ···& 0& 0\\ 0& 1& {\alpha }_2^n& 1& 0& ···& 0& 0\\ ·& & & & & & & ·\\ ·& & & & & & & ·\\ ·& & & & & & & ·\\ 0& ···& & & 0& 1& {\alpha }_{N-1}^n& 1\\ 0& ···& & & & 0& 0& 1\end{array}\right),U_n=\left(\begin{array}{c}\left(\begin{array}{c}{u}_1^n\\ u_2^n\\ ·\\ ·\\ ·\\ u_N^n\end{array}\right)\end{array}\right),V_n=\left(\begin{array}{c}{v}_1^n\\ v_2^n\\ ·\\ ·\\ ·\\ v_N^n\end{array}\right)=\left(\begin{array}{cccccccc}1& 0& 0& & & ···& 0& 0\\ 1& {\alpha }_1^{n-1}& 1& 0& & ···& 0& 0\\ 1& 1& {\alpha }_2^{n-1}& 1& 0& ···& 0& 0\\ ·& & & & & & & ·\\ ·& & & & & & & ·\\ ·& & & & & & & ·\\ 0& ···& & & 0& 1& {\alpha }_{N-1}^{n-1}& 1\\ 0& ···& & & & 0& 0& 1\end{array}\right)\left(\begin{array}{c}{u}_1^{n-1}\\ u_2^{n-1}\\ ·\\ ·\\ ·\\ u_N^{n-1}\end{array}\right)$

其中是A_n第1行的非零元素，是U_n第0行的元素，是V_n第0行的元素，是 A_n-1第1行的非零元素，是U_n-1第0行的元素，其他的以此类推。

把这个三对角矩阵分成不同的块放在每个节点的不同的处理器核上并行，如对一个 12阶的三对角矩阵，可以如下进行划分：

其中a₂,b₂,c₂分别是12阶矩阵的三对角上的第1行的非零元素，其他以此类推。

通过并行分裂法，最终将这个12阶的系数矩阵变成对角形式，如下所示：

其中a₁′～a′₁₂是经过化简后的对角矩阵各行的非零元素。

从而可以并行在不同的处理器上分别求解，从而大大加速了求解过程。

三对角方程的并行加速比如图3所示。通过三级的并行过程，对计算结果与串行程 序进行了比较，得出的计算加速比如图4所示，计算效率如图5所示。可见计算加速比 接近线性，计算效率在32个计算核时还能达到75%以上。

（4）计算的结果进行收集并进行相应整理的步骤。

这个步骤主要是对多级并行计算的结果进行收集并进行相应整理，具体实现过程如 下所示：

①对于上一步骤算出来的计算结果，在收集结果的过程中，要充分考虑计算过程中 结果的形式进行相应的收集。由于采用了交叉步进的计算方法，需要根据交叉步进过程 中二维切片的编号对结果进行整理；

②若对某个二维切片进行了5个交叉步进过程的分解，那么需要按照顺序找到这5 个文件，把结果交叉的放在同一个结果输出文件里，每个文件保存每个二维切片的值；

③在所有二维切片计算完成之后在进行，将所有的二维切片数据组合成三维的数据 形式，最终形成所需要显示的三维数据结果值。

对计算结果进行可视化后的效果图如图6所示。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

以上所述，仅为本发明部分具体实施方式，但本发明的保护范围并不局限于此，任 何熟悉本领域的人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖 在本发明的保护范围之内。