﻿ 一种基于亥姆霍兹分解的大地电磁测深有限元正演预条件解法
 地球物理学报  2019, Vol. 62 Issue (10): 3898-3911 PDF

1. 中国地质大学(北京)地球物理与信息技术学院, 北京 100083;
2. 四川大学灾后重建与管理学院, 成都 610207;
3. 地下信息探测技术与仪器教育部重点实验室, 北京 100083;
4. 深地科学与工程教育部重点实验室, 四川大学, 成都 610065

A Helmholtz decomposition based pre-condition method for magnetotelluric finite element numerical simulation
GUO ZeQiu1,2,4, DONG Hao1,3
1. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China;
2. Institute for Disaster Management and Reconstruction, Sichuan University, Chengdu 610207, China;
3. Key Laboratory of Geo-detection of Ministry of Education, Beijing 100083, China;
4. Key Laboratory of Deep Earth Science and Engineering(Sichuan University), Ministry of Education, Chengdu 610065, China
Abstract: To improve the efficiency of iterative methods solving linear equations in finite element simulation of magnetotelluric sounding,we present a pre-conditioning method exploiting the Helmholtz decomposition of vector fields. Helmholtz decomposition is used to transform the pre-conditioned problem of vector curl-curl equation of time-harmonic Maxwell equation into a Poisson problem based on the vector potential and a Laplacian problem with the scalar potential. With tetrahedral unstructured edge element,the preconditioned problem is discretized on auxiliary node element grid and then solved by algebraic multigrid method (AMG). Then,the COMMEMI synthetic model is used to calculate responses,which is compared with those of integral equation method,and verifies the reliability of the simulation program and the pre-condition method. Additionally,the efficiency of the pre-condition algorithm is tested by discretizing the model with a large range of degrees of freedoms (DoFs). The results show that the new algorithm can significantly improve the convergence of iterative methods that applied to magnetotelluric numerical simulation and presents significantly higher efficiency than the common pre-conditioner with incomplete LU decomposition. Under the circumstances of comparably large DoFs (over 10 million),it manifests apparent advantages over direct solver in efficiency and memory consuming,and also makes the magnetotelluric numerical simulation with a finite element mesh of relatively large size possible on a mobile workstation.
Keywords: Finite element    Helmholtz decomposition    Magnetotellurics    Pre-conditioning    Numerical simulation
0 引言

1 时谐麦克斯韦方程组的变分形式与散度规范问题

 (1a)

 (1b)

 (2a)

 (2b)

 (3)

 (4)

 (5)

n为边界法向方向，ED为边界场值.本文中，我们在求解域的上边界直接给出均匀电场边界条件，而在四周边界和底边界采用大地电磁测深一维正演计算出的电场作为边界条件.借助矢量恒等式，(4)式可转化为下述弱解问题，即寻找EH(curl; Ω)(H(curl)为希尔伯特向量空间中的旋度相容空间)，对于所有vH(curl; Ω)满足变分弱形式：

 (6)

 (7)

 (8)

2 非结构化四面体棱边有限元

 (9)

 (10)

 (11)

 (12)

 (13)

 (14)

 (15a)

 (15b)

3 基于亥姆霍兹分解的预条件矩阵

 (16)

 (17)

 (18)

 (19)

 (20)

 (21)

 (22)

P=A时，P－1Ax=Ix=P－1b，也就是说第一次迭代就获得了准确解.然而，求解上述系统的难度与求解原始方程完全相同，这使得预条件失去了意义.但反过来，如果可以找到一种可以快速近似求解上述系统的方法，迭代算法就可以迅速收敛，获得满足误差要求的解.

 (23)

 (24)

 (25)

 (26)

 (27)

 (28)

 (29)

 (30)

 (31)

 (32)

 图 1 (a) 节点单元空间中的势场向棱边单元空间中的矢量场的梯度运算；(b)节点单元空间中的矢量场分量向棱边单元空间中的投影 Fig. 1 (a) Gradient operation from scalar potential field in nodal space to vector field in edge space; (b) Projection of an orthogonal vector component in nodal space to vector field in edge space

 (33)

 (34)

4 算例 4.1 可靠性验证

 图 2 COMMEMI 3D-2模型(Zhdanov et al., 1997)的电导率分布，(a)中的黑色实线表示正演计算MT传递函数的测线位置 Fig. 2 Conductivity distribution of COMMEMI 3D-2 model (Zhdanov et al., 1997), the solid line in (a) denotes the profile where the MT transfer functions are calculated
 图 3 利用非结构化有限单元离散化的COMMEMI 3D-2模型 (a)网格剖分情况; (b)以q值表示的网格质量. Fig. 3 COMMEMI 3D-2 model discretized with unstructured finite element mesh (a) Scheme of generated mesh; (b) Mesh quality represented by q value.

 图 4 COMMEMI 3D-2模型0.1 Hz的视电阻率及相位响应对比；其中十字为本研究中算法获得的结果，圆圈响应来自积分方程法的结果(Farquharson and Miensopust, 2011) Fig. 4 Comparison of apparent resistivities and phases calculated for the COMMEMI 3D-2 model at a frequency of 0.1 Hz. The crosses are solution obtained using the presented method while the open circles are from Farquharson and Miensopust (2011)
 图 5 COMMEMI 3D-2模型0.001 Hz的视电阻率及相位响应对比；其中十字为本研究中算法获得的结果，圆圈响应来自积分方程法的结果(Farquharson and Miensopust, 2011) Fig. 5 Comparison of apparent resistivities and phases calculated for the COMMEMI 3D-2 model at a frequency of 0.001 Hz. The crosses are solution obtained using the presented method while the open circles are from Farquharson and Miensopust (2011

 图 6 TFQMR算法预条件求解收敛情况对比 其中实线为利用Krylov子空间算法(PCG)求解的亥姆霍兹分解预条件算法，虚线为利用AMG算法求解的亥姆霍兹分解预条件算法，点线为ILU分解预条件算法，注意在此只截取了前1000次迭代的结果. Fig. 6 Convergence comparison of TFQMR method with different pre-conditioners The solid line denotes Helmholtz decomposition pre-conditioning with Krylov space solvers (PCG). Dashed line is the same pre-conditioning with AMG solver, while the dotted line shows result for incomplete LU pre-conditioning. Note that only the residual from first 1000 iterations are shown.

 图 7 不同频率时HD-PCG、HD-AMG与MUMPs的求解效率对比 Fig. 7 Walltime comparison of TFQMR method with different pre-conditioners (PCG, AMG) and MUMPs for different frequencies
4.2 实用性验证

 图 8 COMMEMI 3D-1模型(Zhdanov et al., 1997)的电导率分布，(a)中的虚线表示正演计算MT传递函数的测线位置 Fig. 8 Conductivity distribution of COMMEMI 3D-1 model (Zhdanov et al., 1997), the dotted line in (a) denotes the profile where the MT transfer functions are calculated
 图 9 COMMEMI 3D-1模型0.1 Hz的视电阻率及相位响应对比；其中十字为本研究中算法获得的结果，实线响应来自Ren et al. (2013) Fig. 9 Comparison of apparent resistivities and phases calculated for the COMMEMI 3D-1 model at a frequency of 0.1 Hz. The crosses are solution obtained using the presented method while the solid lines are from Ren et al. (2013)
 图 10 COMMEMI 3D-1模型10 Hz的视电阻率及相位响应对比； 其中十字为本研究中算法获得的结果，实线响应来自Ren et al. (2013) Fig. 10 Comparison of apparent resistivities and phases calculated for the COMMEMI 3D-1 model at a frequency of 10 Hz. The crosses are solution obtained using the presented method while the solid lines are from Ren et al. (2013)

 图 11 COMMEMI 3D-1模型不同自由度模型正演效率对比图示，MUMPs解法因在千万级别自由度时内存需求过大，无法计算，并且其对频率变化不敏感，因而两条曲线基本重合(以虚线表示耗时变化趋势) Fig. 11 Illustration of computing performance for COMMEMI 3D-1 model with different DoFs. For a DoF of 10 million, MUMPs needs so large memory that the result cannot be presented here, and it′s obvious that its efficiency is almost independent of frequency (the trend for walltime is indicated by the dashed line)

5 结论

References