肖调杰,周 峰,郑翾宇,刘 剑,陈 琳,刘 杰,易明宽, 陈旭光,龚春叶,4,杨 博,甘新标,李胜国,左 克

(1.国防科技大学计算机学院,湖南 长沙 410073;2.国防科技大学高端装备数字化软件实验室,湖南 长沙 410073; 3.东华理工大学地球物理与测控技术学院,江西 南昌 330013;4.国家超级计算天津中心,天津 300457)

1 引言

地电磁学中的频率域电磁法是一类重要的地球物理方法,其通过观测由不同频率天然场源或人工场源激发的电磁场在地下介质中产生的感应二次场,并分析其分布规律,进而获得地下电性分布特征及相关电性结构参数[1-4]。频率域电磁法种类繁多,包含大地电磁测深法、可控源音频大地电磁法、海洋可控源电磁法等主动源方法以及大地电磁法、音频大地电磁法、甚低频法等被动源方法,具有探测深度大、分辨率高及抗干扰能力强等优点,在地球深部结构探测、地热勘探、矿产与油气勘查及环境与工程勘探等诸多领域有着广泛的应用[5-10]。

数值方法和线性方程组求解是频率域电磁法数值模拟的2个重要方面。一方面,有限单元法、有限差分法和积分方程法等是当前频率域电磁法中使用最为广泛的3种数值模拟方法[11-17]。其中有限单元法和有限差分法属于微分方程法,基于Maxwell方程的微分形式,求解时需要将全域作为计算区域,会造成较大的未知量和计算量,当模型规模较大时,其计算自由度可迅速增加至数千万、数亿量级甚至更大。积分方程法基于Maxwell方程的积分形式,是数值模拟算法中应用得最早的一种方法。采用积分方程法计算电磁响应时,只需对目标异常体进行剖分,相比于有限单元法和有限差分法,积分方程所需求解的未知量大大减少。由于引入了格林函数,积分方程法避免了微分方程法中场值传递造成的累积误差,并且不需要施加吸收边界条件或截断边界条件,从而避免了截断边界误差,其数值结果具有半解析解的较高精度[18,19]。另一方面,相比于迭代法获得近似解,直接解法可以获得精确解。在三维数值模拟中,线性方程组求解占据了大部分计算量,求解速度直接关系到三维数值模拟的计算速度。目前方程组求解可以分为直接解法和迭代解法。迭代解法通过迭代过程中的不断修正来获得线性方程组的近似解。相比迭代解法,直接解法计算量非常大,且需要大量存储空间,但是对于解存在的线性方程组,直接解法总能够在有限的时间内获得其精确解,且具有非常高的稳定性[20,21]。

当前,频率域电磁积分方程法数值模拟存在内存受限和计算速度慢等问题,而大规模模型多频点、多场点的计算已成为其实际应用中的瓶颈。随着Harrington[22]在1968年提出了矩量法,此后积分方程法在计算电磁领域得到了快速发展。而在三维地电领域,Hohmann[23]在1975年基于六面体网格阐述了奇异核积分技术,开启了将积分方程法用于求解三维地电问题的大门,随后Weidelt[24]推导得到了层状介质张量格林函数的积分表达式,奠定了积分方程法在地电问题中的应用基础。之后,诸多研究人员对积分方程法在频率域电磁法中的应用开展了大量的研究[10,12]。积分方程法在求解过程中需要求解稠密线性方程组,若求解区域的未知量为W,则其存储需量为O(W2),直接求解的计算复杂度为O(W3),随着计算规模的增大,面临着内存受限和计算速度慢等问题。在优化存储方面,Ting等[25]提出利用格林函数的对称性来减少系数矩阵的存储,但是随着模型网格数量的进一步增大,依然面临着内存受限的问题。在提高计算速度方面,Wannamaker等[26]基于大网格单元的格林系数插值出小网格单元的格林系数来加快获得格林系数矩阵;Aksun[27]、袁建生等[28]采用复镜像法加快格林函数的计算。为加快方程组的求解速度,Tripp等[29]把系数矩阵拆分成了若干个小矩阵,Millard等[30]提出了稳定双共轭梯度的快速傅里叶变换法。此外,Zhadnov等[31,32]、Abubakar等[33]、Ueda等[34]和Kruglyakov等[35]深入研究了Born近似、拟线性近似等方法,以降低求解精度为代价提高了计算速度。而随着计算规模、场点和频点的迅速增加,所需内存和计算量将急剧递增。上述研究在一定程度上缓解了内存受限和计算速度慢的问题,但未能从根本上给出合理的解决方案。

因此,本文面向频率域电磁场数值模拟,采用积分方程法、分布式存储和直接解法等技术,设计了一种多层级多粒度混合并行策略,解决了内存受限和计算速度慢等问题,实现了频点间并行、阻抗矩阵并行填充、方程组并行求解的快速、高精度、高可扩展性的频率域电磁三维积分方程数值模拟方法。

2 物理模型

图1所示为频率域电磁法三维地电模型示意图。在地表有一接地导线作为电性发射源(Transmitter),源的长度为L,发射电流为I。积分方程法求解范围是电磁参数区别于背景介质的异常区域,图1中所示背景区域是电导率为σb的均匀半空间,其中存在电导率为σ的异常区域。给定发射源、介质电性参数及空间几何等信息,可通过数值模拟计算得到整个空间的电磁场分布。

Figure 1 Diagram of geo-electromagnetic model图1 三维地电模型示意图

3 积分方程法

研究任何形式的电磁场都必须从Maxwell方程组出发,电磁场通过傅里叶变换,可将时间域的电磁场转换为一系列谐变场的组合,取时谐因子为e-iwt,在准静态条件下,Maxwell方程组变为式(1)～式(4)[5]:

×E=iωμΗ

(1)

×H=(σ-iωε)E+Js

(2)

·εE=ρ

(3)

·μH=0

(4)

其中,E是电场强度矢量,H是磁场强度矢量,ρ是自由电荷密度,ε是介质的介电常数,μ是介质磁导率,σ是介质的电导率,Js是激发源,ω=2πf是角频率,f是频率,i是虚数单位。

对式(1)求旋度后,结合式(2),可得到电场双旋度方程,如式(5)所示:

××E-iωμ(σ-iωε)E=iωμJs

(5)

根据叠加原理,将总电场分解为入射电场Einc和散射电场Esca,有:

E=Einc+Esca

(6)

采用电张量格林函数,并结合式(5),可以得到积分公式,如式(7)所示:

(7)

其中,Ωs是场源积分区域,Ω是异常体积分区域,G是并矢Green函数。当背景模型为均匀半空间时,G是均匀半空间并矢Green函数;当背景模型为层状介质时,G是层状介质并矢函数[25,26]。

在进行数值模拟之前,需要对式(7)进行剖分离散。本文采用四面体单元,将异常体离散成N个四面体子单元。采用适用于介质的体积分方程法,假设每个单元内部的电导率和电场保持不变,取为单元中心ri(i=1,2,…,M)处的值(如图2所示,其中Ex,Ey和Ez分别为x,y和z方向的电场强度),可得式(8):

E(rm)=Einc(rm)+

(8)

其中,Ei为单元中心处的电场。

Figure 2 Diagram of tetrahedral element图2 四面体单元示意图

将式(8)写成矩阵形式,可形成关于异常体区域3N个未知电场的线性方程组,如式(9)～式(12)所示:

A·E=Einc

(9)

(10)

(11)

(12)

求解式(8),可得到异常体区域的电场值。而最终某一点的总场由该点的入射场和散射场叠加而成,再通过式(8)即可获取空间中任一点处的电场值。

4 并行策略

4.1 多级并行架构

在计算之前,采用国产网格剖分软件YH-Grid或开源软件Tetgen对模型进行网格离散,并准备好相应的输入文件。电磁法的频点间计算无数据依赖关系,具有天然的可并行性。对单频点的串行算法进行计算时间占比分析,发现其阻抗矩阵填充、方程组求解和观测点场值计算占总计算时间90%以上。图3为本文设计的并行算法流程示意图。首先,采取进程分组的方式将不同频点划分到不同的进程组。然后,各进程组的主进程读取输入文件并广播给组内所有进程。接着,并行填充阻抗矩阵A,采用分布式存储。最后加入边界条件,对方程组进行并行求解。求解方程组之后,可进一步对观测点的场值进行并行计算。最后进行数据收集和后处理。

Figure 3 Flow diagram of multi-level parallel algorithm图3 多级并行算法流程示意图

4.2 负载均衡

为保持良好的负载均衡与可扩展性,采取分布式存储,采用广泛应用在ScaLAPACK、SLATE、MAGMA等稠密线性代数库中的二维Block- cyclic数据映射方法。假设一共有N个四面体单元,则阻抗矩阵A的大小为3N×3N。将矩阵元素以块的形式不重复地循环分配/存储在P×Q个不同进程之中,其中P和Q分别为二维进程网格中垂直方向和水平方向的进程个数。二维进程网格如图4所示,在逻辑上将9N2个矩阵元素分成一个个大小为mb×nb的分块小矩阵,其中每个格子代表一个分块小矩阵,相同颜色的分块小矩阵存储在同一个进程中,每个进程保存的矩阵元素个数为9N2/(P×Q),一般要求缓存M≤9N2/(P×Q)。

Figure 4 Block-cyclic mapping图4 Block-cyclic映射

设当前进程在进程网格中的索引为(my_row,my_col),图4a中矩阵元素的全局索引为(ig,jg),图4b中当前进程中矩阵元素的局部索引为(il,jl),索引从0开始,则有:

my_row=(ig/mb)%P

(13)

my_col=(jg/nb)%Q

(14)

当前进程局部索引和全局索引之间的关系满足式(15)～式(18):

il=(ig/mb)/P×mb+ig%mb

(15)

jl=(jg/nb)/Q×nb+jg/nb

(16)

ig=P×mb×(il/mb)+il%mb

(17)

jg=Q×nb×(jl/nb)+jl%nb

(18)

4.3 网格分发

根据式(10)可知,阻抗矩阵A由四面体两两之间的相互关系系数矩阵组合而成,如图5所示,每个四面体单元都对应3行和3列。

Figure 5 Impedance matrix A图5 阻抗矩阵A

以总进程为6、P×Q=2×3、阻抗矩阵A分成9×9个分块小矩阵Tij(i=0,…,8;j=0,…,8)为例,每个分块小矩阵大小为NB×NB,Tij如式(19)所示:

(19)

如图6所示,按照块循环(block-cyclic)方式将阻抗矩阵A映射到2×3的进程网格中。

Figure 6 Impedance matrix图6 阻抗矩阵分布示意图

如图7所示,将每行或每列块(block)对应的单元网格信息存放在一个包(bucket)中,矩阵块是方阵,第i行block和第i列block涉及到的网格单元相同,对应同一个bucket,在图中不做展示。其中每个四面体需要保存其全局编号(index)、属性编号(id)、4个节点位置信息(P0(x0,y0,z0),P1(x1,y1,z1),P2(x2,y2,z2),P3(x3,y3,z3)),即2个整数和12个浮点数,分别为(index,id,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3)。

Figure 7 Grid distribution图7 网格分发

根据图6b,进程(0,0)存储的矩阵块如图8所示,进程(0,0)所存储的矩阵块为T0,0,T0,3,T0,6,T2,0,T2,3,T2,6,T4,0,T4,3,T4,6,T6,0,T6,3,T6,6,T8,0,T8,3,T8,6,填充Tij需要第i和第j个bucket,其中行涉及到的bucket为第0,2,4,6和8个bucket,列涉及到的bucket为第0,3和6个bucket。因此,一起需要第0,2,3,4,6和8个bucket。将行和列涉及到的bucket指针分别保存在rowGrid和colGrid2个不同的数组中。一般地,第i个bucket由进程(i%P,i%Q)进行广播。首先,负责读入网格和切分网格到bucket的主进程将第i个bucket发送给进程(i%P,i%Q),再由该进程在row_comm和col_comm分别广播,行进程号等于i%P的进程将收到的bucket指针保存在rowGrid中。列进程号等于i%Q的进程将收到的bucket指针保存在colGrid中,即完成网格分发。

其中,负责广播的进程保存的bucket指针同时保存在rowGrid和colGrid中。本文引入isBoth的哈希表进行判断,以避免内存重复释放。

Figure 8 Matrix block of process (0,0)图8 进程(0,0)存储的矩阵块

4.4 阻抗矩阵填充

网格分发完成之后,各进程已得到所需网格信息。接下来,只需要在rowGrid和colGrid中对每个四面体分别进行计算,然后保存到本地阻抗矩阵即可。假设rowGrid有m个bucket,colGrid有n个bucket,阻抗矩阵填充如算法1所示[36]。

算法1阻抗矩阵填充

输入：m个bucket组成的网格块rowGrid和n个bucket组成的网格块colGrid。

输出：本地阻抗矩阵A。

forj=0:ndo

fori=0:mdo

Aptr=A+i*NB+j*NB*ld;/*跳转到对应NB块的首地址*/

nElement←colGrid的四面体个数;

mElement←rowGrid的四面体个数;

fork=0:nElementdo

计算本地填充的列方向开始位置n0和结束位置n1;

forp=0:mElementdo

计算本地填充的行方向开始位置m0和结束位置m1;

计算四面体之间产生的矩阵T;

Aptr(m0:m1,n0:n1)←

T(m0:m1,n0:n1);

endfor

endfor

endfor

endfor

4.5 复数稠密线性系统LU分解

积分方程法最终形成一个线性稠密的复数线性方程组,根据式(10)可知,采用部分选主元的LU分解进行求解[37],如式(20)所示:

(20)

其中,A11∈CNB×NB,A12∈CNB×(N-NB),A21∈C(N-NB)×NB,A22∈C(N-NB)×(N-NB)。

4.6 观测场值并行计算

已知异常体各离散单元的场值,便可根据式(8)求得任意观测点处的场值。各观测点的场值为背景场与二次场之和,其中二次场为地下异常体在该处产生的感应场值,由异常体中各四面体在该处产生的场值累加组成。当模型规模较大、观测点数较多时,观测点处的场值计算耗时巨大。考虑到实际应用场景,进程数通常少于观测点数,因此在各进程组/单个频点内对观测点处的场值进行并行计算。采用主从模式的并行加速策略,各进程组内主进程负责任务调度和结果收集,从进程负责计算并向主进程返回计算结果。计算流程如图9所示。

Figure 9 Flow chart of parallel calculation of field at observation points图9 观测点场值并行计算流程图

5 验证与测试

5.1 正确性和准确度验证

如图10所示,在电导率为0.02 S/m的均匀半空间中有一个电导率为0.2 S/m的低阻体。低阻体的埋深为100 m,尺寸为120 m × 200 m × 400 m,中心坐标为(0 m,1 000 m,300 m)。偶极子源沿y方向,长度为1 m,发射电流为100 A,频率为3 Hz,中心坐标为(0 m,0 m,0 m)。沿y轴方向在异常体上方布设间距为20 m的41个观测点,起点位置为(0 m,600 m,0 m),终点位置为(0 m,1 400 m,0 m)。

Figure 10 Testing model:There is a conductive body embedded in a half-space图10 测试模型:均匀半空间中存在一低阻体

选取Ey的实部、虚部与前人结果[38]进行对比,如图11所示,可以看到,实部和虚部的相对误差在6%以内。

Figure 11 Comparison of electric field amplitudes图11 电场振幅对比

Figure 13 Distribution of real component of Ey图13 Ey实部分布

5.2 性能测试

以可控源为例,如图12所示。按对数在10-1～102Hz均匀取16个频点,共16 × 12495个未知量。若采用串行算法中的全存储策略而不是本文中的分布式存储,阻抗矩阵所需存储空间大小为16×(12495×2)2×8/10243=74 GB,而随着问题规模的增加,所需存储空间会进一步极速剧增。因此,采用本文算法,在天河高性能计算系统上进行测试,每个节点分配32个进程,共设置861个观测点。测试了1个节点32个进程到256个节点8 192个进程的不同规模。

Figure 12 Large-scale model for testing图12 大规模测试模型

参考Dublin Test Model,设计组合体模型如图12所示。在图12中,电导率σb为0.001 25 S/m的均匀半空间存在一个埋深500 m的低阻异常组合体,组合体的电导率σ1、σ2及σ3分别为0.02 S/m,0.01 S/m和0.005 S/m,对应块体的尺寸分别为300 m×400 m×600 m,300 m×400 m×600 m,2000 m×400 m×400 m,发射场源为沿y方向布设,长度100 m,电流大小为2 A,观测区域位于异常体的正上方,x方向:-2 000 m～2 000 m,y方向:-1 000 m～1 000 m,x和y方向间距100 m,共布设861个观测点。坐标系的源点O设定在异常组合体在水平面投影的中心点处。

Figure 14 Distribution of imaginary component of Ey图14 Ey虚部分布

计算结果如图13和图14所示,分别展示了Ey的实部和虚部,可以看到异常体上方的电场有明显的变化,不同频点对异常体均有所反映。

可扩展性测试结果如表1和图15所示。可以看到,本文开发的并行程序具有较好的可扩展性,相比1个节点32进程,当测试规模达到256节点8 192进程时,加速比为69.69,并行效率为27.22%。

Table 1 Speedup and parallel efficiency of large-scale model表1 大规模模型的加速比及并行效率

Figure 15 Speedup and parallel efficiency of large-scale model图15 大规模模型的加速比和并行效率

6 结束语

本文详细描述了一种频率域电磁积分方程法大规模数值模拟程序的设计、实现方案和测试结果。通过对程序的验证和测试可以看出,该程序具有较高的精度和较好的可扩展性,具备在大型超算平台上大规模运行的能力,可以满足地电磁学中大地电磁、可控源音频大地电磁等精细数值模拟的需求。