一种高效的时域有限元区域分解并行方法与流程

本发明属于大规模并行计算技术领域，特别是一种高效的时域有限元区域分解并行方法。

背景技术：

随着计算机硬件技术的飞速发展，计算规模的增大，单核计算已不能满足目前的需求，多核大规模并行计算成为目前的主流方向。传统的umfpack求解器不能并行求解稀疏矩阵，而且在求解矩阵规模上有一定的限制，而并行mumps求解器虽然可以并行求解，但是随着进程数的增多，效率降得很快，失去了高度并行的特性。

技术实现要素：

本发明的目的在于提供一种高效的时域有限元区域分解并行方法，以提高求解大规模电磁问题的能力。

实现本发明目的的技术解决方案为：一种高效的时域有限元区域分解并行方法，步骤如下：

步骤1，建立求解模型，使用四面体网格对该模型进行离散，得到模型的结构信息，所述结构信息包括四面体的节点信息以及单元信息；

步骤2，消息传递接口mpi初始化，确定进程的总数目，使子区域与进程号一一对应，将单元划分到各个子区域中，单元中的未知量无重复地分配到各个子区域中，通过时间步长大小确定缓冲区厚度，设置仿真参数，建立每个子区域和全局计算区域的编号索引以及相邻子区域之间的通信索引；

步骤3，从以电场强度和磁通量为未知量的一阶麦克斯韦旋度方程出发，对等式两边采用伽辽金法测试，并用基函数展开，得到最终的迭代公式，根据该迭代公式填充各求解子区域的矩阵；

步骤4，通过cn差分格式时域有限元的迭代公式进行时间迭代，一次时间步迭代里有两次通信，分别通信缓冲区内正确的电场值和磁通量，该通信只发生在相邻区域之间，迭代结束得到空间中的电场值和磁通量。

进一步地，步骤2所述消息传递接口mpi初始化，确定进程的总数目，使子区域与进程号一一对应，将单元划分到各个子区域中，单元中的未知量无重复地分配到各个子区域中，通过时间步长大小确定缓冲区厚度，设置仿真参数，建立每个子区域和全局计算区域的编号索引以及相邻子区域之间的通信索引，具体为：

根据进程总数目将时域有限元计算区域划分为与进程数相同个数的子区域，每个单元都有对应所在的进程号，而相邻子区域交界处的未知量信息则要无重复地分配到各个子区域中，不能存在两个相邻子区域共用未知量的情况，对每个子区域根据时间步长的大小进行区域扩展形成缓冲区，对每个子区域的缓冲区的电场和磁场进行整合编码，建立相邻子区域之间的通信索引。

进一步地，步骤3所述从以电场强度和磁通量为未知量的一阶麦克斯韦旋度方程出发，对等式两边采用伽辽金法测试，并用基函数展开，得到最终的迭代公式，根据该迭代公式填充各求解子区域的矩阵，具体如下：

基于电场强度和磁通量的一阶麦克斯韦旋度方程如式(1)、(2)：

式中，t代表时间，是拉普拉斯算子，ε、μ分别表示离散单元的介电常数和磁导率，和分别表示电场强度和磁通量；

对式(1)采用whitney-i型的棱边基函数进行伽辽金法测试，对式(2)采用whitney-ii型的面基函数进行伽辽金法测试，并将和分别用whitney-i型棱基函数和whitney-ii型面基函数展开，最终得到：

其中，e和b分别代表求解区域的电场值和磁通量，[t^e]、[p^e]、[k^e]、[q^e]分别为形成的体积分矩阵，各矩阵块如下：

因为基函数在四面体单元内的线性组合关系，体现在[t^e]、[k^e]、[p^e]的转化关系上：

[p^e]＝[c^e]^t[k^e]^t＝[c^e]^t[k^e](5)

通过topological转换矩阵将式(5)带入式(3)、(4)得到：

进一步地，步骤4所述通过cn差分格式时域有限元的迭代公式进行时间迭代，一次时间步迭代里有两次通信，分别通信缓冲区内正确的电场值和磁通量，该通信只发生在相邻区域之间，迭代结束得到空间中的电场值和磁通量，具体如下：

将式(6)和(7)分别按照cn差分格式展开得到：

[t^e]{e}ⁿ⁺¹＝[t^e]{e}ⁿ+δt[c^e]^t[k^e]{b}ⁿ(8)

{b}ⁿ⁺¹＝{b}ⁿ-δt[c^e]{e}ⁿ⁺¹(9)

其中，[t^e]、[c^e]、[k^e]分别为形成的矩阵；

求出各个子区域的场值之后，各个子区域之间相互通信，通信内容为缓冲区的正确电场值和磁通量；

整个计算区域内所有的单元在每个时间步，均根据所属的子区域进行电场值和磁通量更新，空间中每点均得到正确的场值，至此一个时间步计算全部完成；重复上述步骤直至时间迭代结束。

本发明与现有技术相比，其显著优点为：(1)在传统时域有限元法基础上引入了因果区域分解方法，在通过增加子区域的个数大大提高求解问题规模的同时又可以高效地并行计算，实现隐格式时域有限元法的大规模并行计算；(2)有效地节省时域有限元法仿真计算时间，提高其求解大规模电磁问题的能力，在实际应用中具有重大意义。

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

附图说明

图1为计算区域划分示意图。

图2为本发明实施例中的金属谐振腔结构示意图。

图3为本发明实施例中金属谐振腔的时域波形比较图。

具体实施方式

本发明高效的时域有限元区域分解并行方法，步骤如下：

步骤1，建立求解模型，使用四面体网格对该模型进行离散，得到模型的结构信息，所述结构信息包括四面体的节点信息以及单元信息；

步骤2，mpi(消息传递接口)初始化，确定进程的总数目，对每个进程进行编号，使子区域与进程号一一对应，将四面体单元按照中心点坐标划分到不同的区域中，通过时间步长大小确定缓冲区厚度，设置仿真参数，执行程序前处理；建立每个子区域和全局计算区域的编号索引以及相邻子区域之间的通信索引。

根据进程总数目将时域有限元计算区域划分为与进程数相同个数的子区域，每个单元都有对应所在的进程号，而相邻子区域交界处的未知量信息则要无重复地分配到各个子区域中，不能存在两个相邻子区域共用未知量的情况，对每个子区域根据时间步长的大小进行区域扩展形成缓冲区，对每个子区域的缓冲区的电场和磁场进行整合编码，建立相邻子区域之间的通信索引。

图1中的实线为步骤1中的计算区域划分的边界，每个区域向外扩充一部分缓冲区，即图1中阴影部分区域，区域1的实际计算区域为位于计算区域2、3中虚线所共同围成的计算区域。对每个区域实线与虚线所围成的缓冲区域内的未知量进行编码，建立相邻计算区域之间的通信索引；

步骤3，从以电场强度和磁通量为未知量的一阶麦克斯韦旋度方程出发，对等式两边采用伽辽金法测试，并用基函数展开，得到最终的迭代公式，根据该迭代公式填充各求解子区域的矩阵；

基于电场强度和磁通量的一阶麦克斯韦旋度方程如式(1)、(2)：

式中，t代表时间，是拉普拉斯算子，ε、μ分别表示离散单元的介电常数和磁导率，和分别表示电场强度和磁通量。

对式(1)采用whitney-i型的棱边基函数进行伽辽金法测试，对式(2)采用whitney-ii型的面基函数进行伽辽金法测试，并将和分别用whitney-i型棱基函数和whitney-ii型面基函数展开，最终得到：

其中，e和b分别代表求解区域的电场值和磁通量，[t^e]、[p^e]、[k^e]、[q^e]分别为形成的体积分矩阵，各矩阵块如下：

其中[t^e]为whitney-i型的棱边基函数形成的质量矩阵，[k^e]为whitney-ii型的面基函数形成的质量矩阵。对whitney-i型的棱边基函数取旋度，就会发现它们是whitney-ii型的面基函数的线性组合，即：在四面体单元内这种线性组合关系体现在[t^e]、[k^e]、[p^e]的转化关系上：

[p^e]＝[c^e]^t[k^e]^t＝[c^e]^t[k^e](5)

其中，是topological转换矩阵，仅仅含有元素0，1，-1，它是一个离散的旋度操作，体现的是对whitney-i型的棱边基函数做偏微分离散与whitney-ii型的面基函数的映射关系。将式(5)带入式(3)、(4)得到：

步骤4，通过cn差分格式时域有限元的迭代公式进行时间迭代，一次时间步迭代里有两次通信，分别通信缓冲区内正确的电场值和磁通量，该通信只发生在相邻区域之间，迭代结束得到空间中的电场值和磁通量。

将式(6)和式(7)分别按照隐式(cn)差分格式展开得到：

[t^e]{e}ⁿ⁺¹＝[t^e]{e}ⁿ+δt[c^e]^t[k^e]{b}ⁿ(8)

{b}ⁿ⁺¹＝{b}ⁿ-δt[c^e]{e}ⁿ⁺¹(9)

其中，[t^e]、[c^e]、[k^e]分别为形成的矩阵；

将求解区域的未知量按照上述迭代公式进行迭代之后，每次迭代分别通信缓冲区正确的电场值和磁通量，空间中每点均可得到正确电场值，至此一个时间步中电场的计算全部完成。每个时间步重复上述步骤，即可得到所需的结果。最后释放空间，mpi并行结束多个进程。

实施例1

为了验证本发明的正确性与有效性，下面分析了一个金属谐振腔。

一个如图2所示的金属谐振腔，尺寸为1m*1m*1m，在z方向0.2m处加源，加源电场方向为y方向，观察频率0-0.4ghz，z方向0.8m处设置观察点。采用0.05m四面体剖分，剖出来总共约3万个体，未知量约18万。时间步为1000步，时间步长为8.33ps，缓冲区每个进程向外扩0.075m的距离。图3为时域有限元区域分解并行方法采用16个进程计算出来的时域波形和传统时域有限元采用并行mumps求解器计算出来的时域波形比较，验证了本方法的正确性。表1给出了本专利提出方法与时域有限元法采用传统的并行mumps求解器时计算时间对比，从表中可以看出本发明大大减少了计算时间。从并行效率中可以看出本发明实现了时域有限元区域分解并行方法的高度并行。

表1.计算时间对比

本发明提到的方法求逆采用的是umfpack求解器，在求解不同规模的矩阵时所消耗的时间不是线性关系的。以两个进程为例，单进程求解一个大矩阵和两个进程分别求解两个小矩阵的时间大于两倍关系，而所分进程数较少，进程间通信的未知量相对进程内求解的总未知量较少时，通信所消耗的时间对总的求解时间影响不大，所以出现了并行效率超过100％的情况。