A hybrid boundary element-finite element approach to modeling plane wave 3D electromagnetic induction responses in the Earth

Journal of Computational Physics 258 (2014) 705–717

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A hybrid boundary element-finite element approach to modeling plane wave 3D electromagnetic induction responses in the Earth Zhengyong Ren a,b,∗ , Thomas Kalscheuer a , Stewart Greenhalgh a , Hansruedi Maurer a a b

Institute of Geophysics, Department of Earth Sciences, ETH Zurich, Zurich, 8092, Switzerland School of Geosciences and Info-physics, Central South University, Changsha, 410083, Hunan, China

a r t i c l e

i n f o

Article history: Received 27 November 2012 Received in revised form 31 October 2013 Accepted 6 November 2013 Available online 12 November 2013 Keywords: Boundary element-finite element Adaptive multi-level fast multipole Magnetotelluric Electromagnetic induction Radio-magnetotelluric Arbitrary topography

a b s t r a c t A novel hybrid boundary element-finite element scheme which is accelerated by an adaptive multi-level fast multipole algorithm is presented to simulate 3D plane wave electromagnetic induction responses in the Earth. The remarkable advantages of this novel scheme are the complete removal of the volume discretization of the air space and the capability of simulating large-scale complicated geo-electromagnetic induction problems. To achieve this goal, first the Galerkin edge-based finite-element method (FEM) using unstructured meshes is adopted to solve the electric field differential equation in the heterogeneous Earth, where arbitrary distributions of conductivity, magnetic permeability and dielectric permittivity are allowed for. Second, the point collocation boundary-element method (BEM) is used to solve a surface integral formula in terms of the reduced electrical vector potential on the arbitrarily shaped air–Earth interface. Third, to avoid explicit storage of the system matrix arising from large-scale problems and to reduce the horrendous time complexity of the product of the system matrix with an initial vector of unknowns, the adaptive multilevel fast multipole method is applied. This leads to a matrix-free form suitable for the application of iterative solvers. Furthermore, a highly sparse problemdependent preconditioner is developed to significantly reduce the number of iterations used by the iterative solvers. The efficacy of the presented hybrid scheme is verified on two synthetic examples against different numerical techniques such as goal-oriented adaptive finite-element methods. Numerical experiments show that at low frequencies, where the quasi-static approximation is applicable, standard FEM methods prove to be superior to our hybrid BEM–FEM solutions in terms of computational time, because the FEM method requires only a coarse discretization of the air domain and offers an advantageous sparsity of the system matrix. At radio-magnetotelluric frequencies of a few hundred kHz, the hybrid BEM–FEM scheme outperforms the FEM method, because it avoids explicit storage of the system matrices as well as dense volume discretization of the air domain required by FEM methods at high frequencies. In summary, to the best of our knowledge, this study is the first attempt at completely removing the air space for large scale complicated electromagnetic induction modeling in the Earth. © 2013 Elsevier Inc. All rights reserved.

*

Corresponding author at: School of Geosciences and Info-physics, Central South University, Changsha, 410083, Hunan, China. E-mail addresses: [email protected] (Z. Ren), [email protected] (T. Kalscheuer), [email protected] (S. Greenhalgh), [email protected] (H. Maurer). 0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.11.004

706

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

1. Introduction Plane wave electromagnetic induction methods have become increasingly popular to study the lithospheric structure of the Earth [1] or to investigate the shallow subsurface in engineering and environmental applications [2]. To accurately invert field data for complicated distributions of the electrical parameters of the Earth, there is the need for accurate and efficient algorithms to simulate plane wave electromagnetic induction responses. Currently, there are three numerical techniques which can simulate electromagnetic induction within a complicated 3D heterogeneous Earth, namely finite-difference methods (FDM) [3,4], finite-element methods (FEM) [5–8] and volume integral methods [9,10]. Unfortunately, finite-difference and finite-element methods need to include the volume discretization of the air space in simulations [11]. There are two benefits from removing the volume discretization of the air space from the computational domain. Firstly, the computational domain can be reduced to be predominantly the heterogeneous Earth, which would dramatically decrease the number of unknowns. Secondly, the potential numerical errors arising from the volume discretization of the air space can be eliminated. At high frequencies, the second benefit becomes more critical because here displacement currents play an important role [12]. The propagation characteristics of the electromagnetic fields become so critical that rather dense grids in the air space are required to guarantee acceptable accuracy. The reason is that the wavelengths of the electromagnetic fields in the air space are rather small. To adequately approximate the significant variations of the waves in the air space, careful mesh design in the air space is required. Therefore, the removal of the volume discretization of the air space will avoid serious numerical errors due to improper mesh discretization strategies. In our previous work [13], a boundary element approach was successfully developed for handling arbitrary surface topography by utilizing surface discretization of the air–Earth interface. However, this BEM cannot deal with the case of arbitrary distributions of the electrical parameters in the subsurface. Therefore, we introduced a finite-element modeling code [14] that deals with arbitrarily complicated models. Combining the benefits of the above two methods, we present here a hybrid BEM–FEM scheme. This hybrid scheme substitutes the volume discretization in the air space by the surface discretization of the air–Earth interface having arbitrary shape and allows for complicated distributions of electrical parameters in the Earth through volume discretization of the subsurface with a FEM. Unlike other hybrid schemes which were mainly developed for the quasi-static approximation [15,16] or purely high frequency wave propagation problems [17], our hybrid scheme aims to work for a broadband frequency range such as magnetotelluric [18] and radio-magnetotelluric problems [19]. We apply the edge-based finite-element method to solve the electric field equation in the heterogeneous Earth which is discretized by unstructured meshes. Based on the homogeneous electrical properties of the air space, we introduce a surface integral formula which treats the reduced electrical vector potential as its primary variable on the air–Earth interface. The simple but efficient point collocation boundary-element method is used to solve this surface integral formula. These finite-element and surface integral systems are strongly coupled by continuity conditions of the tangential components of both the electric and magnetic fields on the air–Earth interface. To extend its capability to deal with large-scale practical cases [17], an accelerated fast multipole method (FMM) [20] which was specially designed for nonuniform unstructured grids [13] is adopted. Based on the methodology of panel clustering [21], this FMM algorithm can dramatically reduce both the memory cost and computation time on unstructured grids. The core role of this acceleration algorithm is to efficiently compute the product of the system matrix with an initial vector of unknowns, which can be naturally embedded in an iterative solver and leads to a matrix-free scheme that avoids explicit storage of the system matrix. In addition, to accelerate the convergence rate of iterative solvers, a highly sparse and efficient preconditioner is developed. In the following, we briefly present the theory and two numerical examples to verify and demonstrate our hybrid scheme. 2. Methods 2.1. Hybrid scheme Three-dimensional (3D) geo-electromagnetic induction modeling requires specification of the geometrical model configuration as shown in Fig. 1. The domain is composed of the air space Ω0 and the inhomogeneous Earth domain Ω1 which are connected by the air–Earth interface Γ . A plane wave (Einc , Hinc ) is incident on the outer surface enclosing the domains Ω0 and Ω1 , ∂(Ω0 + Ω1 ) = Γ0 + Γ1 . The electrical parameters for Ω0 and Ω1 are σ0 , ε0 , μ0 and σ1 , ε1 , μ1 , respectively, with σ1  σ0 . In a heterogeneous Earth, the parameters are functions of position. We take the positive z-direction as downwards. Our purpose is to compute the distribution of electromagnetic fields in Ω0 and Ω1 , in terms of the known electromagnetic boundary conditions on surfaces Γ0 and Γ1 . The Maxwell equations in the frequency domain (with a time varying factor exp−i ωt ) in the heterogeneous Earth Ω1 can be written as [22]

∇ × E1 = −ξ1 H1 ,

(1)

∇ × H 1 = χ1 E 1 ,

(2)

∇ · χ1 E 1 = 0 ,

(3)

∇ · ξ1 H 1 = 0 ,

(4)

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

707

Fig. 1. Geometry of the 3D electromagnetic induction problem in the global coordinate system (x, y , z) with an arbitrary incident plane wave (Einc , Hinc ). ˆ vˆ and n. ˆ The outward unit normal vector A local coordinate system on the air–Earth interface Γ is introduced, with three unit direction vectors being u, ˆ 0 and nˆ 1 on Γ1 . The volume discretization of the dashed air space Ω0 is removed from the computational domain. The computational domain on Γ0 is n only contains the heterogeneous Earth Ω1 and the surfaces Γ , Γ0 and Γ1 . Boundary conditions for analytical magnetic fields of layered Earth models [28] are applied on the surface ∂(Ω0 + Ω1 ) = Γ0 + Γ1 .

where E1 is the electric field strength (V/m), H1 is the magnetic field (A/m), ξ1 is the impedivity, ξ1 = −i ωμ1 , χ1 is the √ admittivity, χ1 = σ1 − i ω 1 , i = −1 is the imaginary unit, μ1 is the magnetic permeability (H/m), ω = 2π f is the angular frequency, f is the frequency, σ1 = 1/ρ1 is the conductivity distribution of the heterogeneous Earth in S/m (where ρ1 is resistivity in m) and 1 is the dielectric permittivity (F/m). Taking the curl of Eq. (1) and eliminating the magnetic field, we obtain a curl–curl differential equation for the electric field

∇×

1

ξ1

∇ × E1 + χ1 E1 = 0.

(5)

The solutions E1 in the heterogeneous Earth which belong to the Hilbert space with finite energy, H 1 (curl, Ω1 ){V, where |V|2 < inf, |∇ × V|2 < inf}, satisfy the weak form of solution of the partial differential equation. This can be stated that to find E1 ∈ H 1 (curl, Ω1 ) such that for arbitrary Galerkin test functions V ∈ H 1 (curl, Ω1 )





V· ∇ ×

1

ξ1

 ∇ × E1 + χ1 E1 dv = 0.

(6)

Ω1

Using the first vector Green’s theorem [22] on Eq. (6), we obtain the following weak form of solution





1

ξ1

 ∇ × V · ∇ × E1 − k21 V · E1 dv =

ˆ × H1 ) ds, V · (n

(7)

∂Ω1 =Γ +Γ1

Ω1

ˆ is the outgoing unit vector on surfaces enclosing the heterogeneous Earth and k1 is the wavenumber, k21 = −ξ1 χ1 . where n In finite-element methods, the weak form of solution given by Eq. (7) is applied in both the air and Earth domains. At very low frequency, the wavenumber of the air space approaches zero. This causes the weak form of solution of the electric field in the air domain to become unstable, and the uniqueness of the solution cannot be guaranteed [23]. Physically, the reason for this instability is that the explicit condition of vanishing divergence of the current density is ignored in the curl–curl electric field equation (5). In the hybrid scheme, the curl–curl partial differential equation is applied in the heterogeneous Earth which has a finite conductivity such that small values of the wavenumber at low frequency and the associated numerical instabilities can be avoided. Furthermore, we can observe that the coefficients in front of the integrand terms in Eq. (7) are quite different to each other in terms of magnitudes, which could potentially lead to an ill conditioning of the final system of linear equations. To avoid this potential problem, we define a new variable, the reduced version of the electric vector potential T [5], to stabilize the weak form of solution in Eq. (7) T∗0,1 = i ωμ0,1 H = −ξ0,1 H,

ξ1 = μr ξ0 (ξ0 = −i ωμ0 is constant).

Then, the weak form of solution given by Eq. (7) becomes



1 

μr Ω1

 ∇ × V · ∇ × E1 − k21 V · E1 dv −

 Γ





ˆ × T∗1 ds = − V· n

 Γ1

(8)





ˆ 1 × T∗1 ds, V· n

(9)

708

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

ˆ is the outward unit normal vector on the surface Γ (see Fig. 1). where 1  μr  1.1 in the heterogeneous Earth [24], and n Due to the homogeneous nature of the air space, surface integral formulae can be introduced to avoid the volume discretization of the air space. This does not only avoid the very dense meshes at high frequency in the air space, but also avoids numerical errors arising from the volume discretization of the air space. In addition, the dimension of the problem in the air space is reduced by one so that the total number of unknowns is dramatically decreased. In the air space Ω0 , the Maxwell’s equations can be written as

∇ × E0 = −ξ0 H0 ,

(10)

∇ × H 0 = χ0 E 0 ,

(11)

∇ · E 0 = 0,

(12)

∇ · H0 = 0.

(13)

At low frequency, the quasi-static assumption in the air space has two implications. Firstly, the vanishing divergence of both the electric field and the magnetic field means that the vector curl–curl equation for the electric field or the magnetic field can be reduced to three scalar Laplace equations in which the electric field and magnetic field are fully decoupled [15]. Secondly, since the curl of the magnetic field is zero (i.e. χ0 = 0 at low frequency), a simple magnetic scalar potential which satisfies the Laplace equation can be used to represent the magnetic field in the air space [25]. Due to the stability of the Laplace operator, robust hybrid boundary-element finite-element schemes based on potentials can be introduced to solve the low frequency instability of geo-electromagnetic problems [26]. However, if displacement currents must be considered at high frequency, the above simplifications arising from the quasi-static assumption are not valid anymore. To resolve this problem, combinations of scalar and vector potentials were introduced for hybrid boundary-element finite-element schemes. Based on the vanishing divergence of the magnetic flux density (∇ · μH = 0), the A − φ approach utilizes a magnetic vector potential and an electric scalar potential in both the air space and the Earth domain [16]. Similarly owing to the vanishing divergence of the total current density (conduction currents plus displacement currents), the T − Ω approach with a combination of a magnetic vector potential and an electric scalar potential was introduced [27]. However, in these hybrid schemes, the primary variables involve several differential operators so that not only higher order shape functions are required to retrieve the electric and magnetic fields, but also elaborate efforts are needed to enforce the continuity conditions on the air–Earth interface [27]. In this study, to have a generalized hybrid scheme, which works for both low and high frequency and has the benefit of easy enforcement of the continuity conditions on the air–Earth interface, we treat the magnetic field in the air space as the primary variable. Using the vector Green’s theorem [22], we obtain the surface magnetic field integral equation (MFIE, [17]) on the smooth air–Earth interface Γ

1 2





 (nˆ × H0 ) × ∇ G + nˆ × (∇ × H0 )G + (nˆ · H0 )∇ G ds,

H0 (´r) = −

(14)

r∈∂Ω0

ˆ is the outward unit normal vector on ∂Ω0 , the derivative operator is defined for r ∈ ∂Ω0 , r´ is where ∂Ω0 = Γ + Γ0 , n an arbitrary observation point on Γ , r is the source point, and G is the full space Green’s functions of the 3D Helmholtz equation in the air space, i.e. G=

e ik0 R 4π R

,

R = |r − r´ |.

(15)

Using the definition of the reduced electrical vector potential T∗ given by Eq. (8), the above MFIE equation becomes

1 ∗ T (´r) = − 2 0













ˆ × T∗0 × ∇ G + k20 nˆ × (E0 )G + nˆ · T∗0 ∇ G ds + T∗Γ0 (´r), n

(16)

r∈Γ

ˆ is the unit normal vector defined on the air–Earth interface (see Fig. 1). In the above equation, T∗Γ0 (´r) is the known where n reduced electric vector potential at the observation point r´ ∈ Γ contributed by the known electromagnetic fields over the upper part Γ0 of the outer surface. On the upper surface Γ0 , the primary electromagnetic fields are given as the boundary conditions, which are usually the analytical solutions of 1D layered-Earth models [28]. The assumption of these boundary conditions is based on the fact that if Γ0 is far enough away (such as 5 to 10 wavelengths) from the anomalous bodies, the secondary fields on Γ0 can be ignored with acceptable accuracy [6,7]. On the air–Earth interface Γ , the boundary conditions for the electric field and the reduced vector electric potential are

ˆ × E0 = nˆ × E1 , n ∗

∗

ˆ × T0 = nˆ × T1 . n

(17) (18)

Based on Huygen’s principle [29], the tangential component of the reduced vector electric potential is referred to as the ˆ × T∗ and its normal component is referred to as the virtual electric charge density virtual electric current density J = n

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

709

β = nˆ · T∗ . Referring to Eqs. (9) and (16), it means that the electric field in Ω1 is physically generated by the virtual surface electric currents and the virtual surface electric currents are related to the virtual surface charges. At low frequencies, i.e. where k20 → 0, Eq. (16) states that on Γ the virtual electric surface charges and currents are independent of the electric field under the ground. The small variations of the electric field due to anomalous bodies could not be observed at low frequencies. To avoid this decoupled behavior, we replace the surface electric charge by the surface electric field distribution, that is

β = −∇ · (nˆ × E).

(19)

ˆ (´r)× to Eq. (16), we obtain a strongly coupled hybrid scheme (named E-J Then applying the cross product operation n scheme) between the electric fields and the reduced vector electric potentials (magnetic fields) as



1 

 ∇ × V · ∇ × E − k21 V · E dv −

μr Ω1

1 2



ˆ (´r) × J(´r) + n



 V · J ds = −

V · JΓ1 ds,

(20)

ˆ × EG − ∇ · (nˆ × E)∇ G ds = JΓ0 (´r). J × ∇ G + k20 n

(21)

Γ



Γ1



r∈Γ

2.2. FEM and BEM discretizations

N

The heterogeneous subsurface Ω1 is discretized into disjoint source tetrahedra, Ω1 = s=t 1 Ω1s , where Ω1s is the sth source tetrahedron, N t is the total number of tetrahedra. The air–Earth interface Γ is discretized into N tΓ source triangles,

Γ =

NtΓ

s=1 s . Here, source tetrahedra and triangles are referred to as elements in which fields are approximated by shape functions. Observation elements (which will be used later) are the elements in which test functions are defined. The virtual surface current density J is represented in a local coordinate frame (u , v , n) on Γ (see Fig. 1), that is

ˆ + J v vˆ , J = Juu

ˆ = uˆ × vˆ . n

(22)

Therefore, the surface formula in the E-J scheme can be decomposed into two scalar equations

1 2 1 2











ˆ × EG − ∇ · (nˆ × E)∇ G ds = uˆ (´r) · JΓ0 (´r), J × ∇ G + k20 n

ˆ (´r) · nˆ (´r) × J u (´r) + u r∈Γ



ˆ × EG − ∇ · (nˆ × E)∇ G ds = vˆ (´r) · JΓ0 (´r). J × ∇ G + k20 n

ˆ (´r) × J v (´r) + vˆ (´r) · n

(23)

(24)

r∈Γ

The electric field in the Earth is expanded as

E=

Ne

E s Ns .

(25)

s =1

Here, Ns is the divergence-free first order Nedelec shape function [30] on the sth edge, E s is the circulation of electric fields on the sth edge, N e is the total number of edges in the domain Ω1 . On the air–Earth interface, we assume that in each triangle, the virtual surface current density J is constant. It means that J u and J v can be approximated in terms of zero-order shape functions on Γ , that is Γ

Ju =

Nt

J us δs ,

(26)

J sv δs ,

(27)

s =1 Γ

Jv =

Nt

s =1

where N tΓ is the total number of triangles on the air–Earth interface Γ , and δ is the Dirac delta function. Substituting Eqs. (25), (26) and (27) into Eqs. (20), (23) and (24), and by integrating the residuals over each observation tetrahedron in Ω1 and each observation triangle on Γ , we obtain a system of linear equations

 AX =

A0 B0 C0

A1 B1 C1

A2 B2 C2

 ¯  E F0 ¯Ju = F1 = F, ¯J v F2

(28)

710

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

where A is the system matrix and X is the vector of unknowns consisting of E¯ = { E s , s = 1, 2, . . . , N e }, ¯Ju = { J us , s = 1, 2, . . . , N tΓ } and ¯J v = { J vs , s = 1, 2, . . . , N tΓ }. The size of the above system of linear equations is ( N e + 2N tΓ ) × ( N e + 2N tΓ ), A0 is a square sparse matrix from the finite-element method part, A1 and A2 (which arise from the FEM to BEM part) are partially sparse matrices of size N e × N tΓ with less than two nonzero entries per row. For a good quality unstructured grid, the maximum number of nonzero entries per row of AFEM = [A0 , A1 , A2 ] is less than 50 [14]. The BEM–FEM matrices B0 and C0 are partially sparse matrices, with a maximum number of nonzeros entries per row of about 1.5 N tΓ . The BEM–BEM matrices B1 , B2 , C1 and C2 are dense, each of size N tΓ × N tΓ . The explicit forms for elements in the above system of linear equations are



1 

A os 0 =

 ∇ × No · ∇ × Ns − k21 No · Ns dv ,

μr

(29)

Ω1



A os 1

=−

ˆ δs ds, No · u

(30)

No · vˆ δs ds,

(31)

Γ

 A os 2 =− Γ

 



ˆ· B os 0 = −δo v Γ

B os 1









ˆ × Ns G − ∇ · nˆ × Ns ∇ G ds , k20 n 

1

= δo δs − δo vˆ ·

(δs uˆ ) × ∇ G ds,

2

(33)

Γ

 ˆ· B os 2 = −δo v

(δs vˆ ) × ∇ G ds,

(34)

Γ

  ˆ· C 0os = δo u

(32)











ˆ × Ns G − ∇ · nˆ × Ns ∇ G ds , k20 n

(35)

Γ

 ˆ· C 1os = δo u

(δs uˆ ) × ∇ G ds, Γ

C 2os =

1 2

(36)



δo δs + δo uˆ ·

(δs vˆ ) × ∇ G ds.

(37)

Γ

The expressions for source terms are

 F 0s = −

Ns · JΓ1 ds,

(38)

Γ1

ˆ · JΓ0 , F 1o = δo u

(39)

= δo vˆ · JΓ0 .

(40)

F 2o

In the above, the superscript o stands for shape functions defined in observation tetrahedra or triangles and the superscript s is for shape functions defined in source tetrahedra or triangles. The volume integrals in A0 are evaluated using analytical forms [22]. Normal Gauss quadrature rules are used to accurately evaluate the regular integrals in A1 and A2 . When the observation triangle overlaps the source triangle, then the surface integrals in B0 , B1 , B2 , C0 , C1 and C2 have weak and strong singularities. A singularity extraction technique is used to accurately evaluate singular integrals [31,13]. 2.3. Evaluation of AX Due to the dense structures of the BEM–BEM matrices B1 , C1 , B2 and C2 in Eq. (28), direct solvers require at least a time complexity of O (( N tΓ )3 ) and a memory cost of O (( N tΓ )2 ). By contrast, iterative solvers require at least a time complexity of O (( N tΓ )2 ) to solve the problem. When the number N tΓ of triangles on the air–Earth interface is increasing, these time complexities cause a bottleneck for large-scale problems. To compute the matrix–vector product AX in Eq. (28), we need to evaluate the following core integral over the air–Earth interface Γ ,

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

711

 J × ∇ G (r, r´ ) ds,

I(´r) =

(41)

Γ

where r´ are the center points of observation triangles, and r are the center points of source triangles on Γ which are discretized into a set of triangles. The adaptive multilevel fast multipole (AMFM) method can reduce the above time complexity of evaluating AX to O ( N tΓ log N tΓ ) accompanied by a linear complexity of memory cost O ( N e + N tΓ ) [13]. Utilizing an adaptive multi-level tree construction algorithm [13] and the panel-clustering idea [21], triangles on the air–Earth interface Γ are grouped into clusters of different levels on a cluster tree. Each cell or element of the tree is a cluster. Each cluster might have a parent cluster at the next higher level on the tree or have several child clusters at the next lower level. On this cluster tree, there are no empty clusters which contain zero triangles. Clusters having no child cluster are referred to as leaf clusters. The core integrals I(´r) are evaluated on the leaf clusters. For each leaf cluster, its integral triangles on the air–Earth interface Γ can be further grouped into near-field integral clusters and far-field integral clusters. The integrals over the near-field clusters are evaluated directly in terms of analytical singularity extraction techniques [31,13]. The far-field clusters can be recognized as a sub-tree structure from the whole cluster tree. The integrals over the far-field clusters then are evaluated from the parent cluster of these far-field clusters, that is

I(´r) = I(´r)near + I(´r)far . Taking advantage of the above adaptively constructed cluster tree structure, the matrix–vector product AX is efficiently computed. In matrix notation, given the initial electric field circulations on edges of Ω1 and surface currents on Γ , the multiplication of the system matrix with the initial field X0 can be written as

AX0 =



AFEM far

Anear BEM + ABEM

X0 ,

(42) far

where AFEM = [A0 , A1 , A2 ], Anear r)near and ABEM X0 is BEM X0 is a matrix–vector product notation containing the core integral I(´ computed through the core integral I(´r)far . The source terms F in Eq. (28) are computed by the standard FMM algorithm [32]. 2.4. An efficient preconditioner The AMFM algorithm is used to accelerate the multiplication of the system matrix with an initial vector as typically employed in iterative solvers [33,34] such as the bi-conjugate gradient stabilized method (BICGSTAB). Due to a combination of two different systems, the condition number of the hybrid scheme might be high. The convergence rate of iterative solvers is quite critical, because it is linearly related to the performance of the AMFM method. Therefore, we need to find an efficient preconditioner to reduce the number of iterations. The ideal preconditioner of a regular matrix A is P = A−1 , so that the preconditioned system would be convergent in one iteration. However, this preconditioner is impractical because it is equal to solving the original system using direct solvers. The difficulty of designing a cheap but efficient preconditioner is arising from the fact that the submatrices B1, B2, C1 and C2 in Eq. (28) are dense. There are several preconditioners available for dense matrices arising from BEM methods, such as the popular incomplete LU preconditioner and block diagonal preconditioners [35–37]. Based on the fact that the full space Green’s function decreases with increasing distance, the combination of the original sparse matrix AFEM from the FEM part and the sparse matrix Anear BEM from the near field integral triangles should be a good choice of preconditioner, i.e.

Pnear =

AFEM Anear BEM



(43)

.

However, because nonuniform grids are used in this study, the local mesh densities might be rapidly varying. This implies near that the number of nonzero entries per row in Anear is not equally distributed and the sparsity of the preconditioner BEM ∈ P is decreased. To avoid this undesired property, a fixed maximum number m of near field triangles is given to generate a well-behaved banded sparse matrix from Anear BEM . Therefore, the preconditioner used in this study is

⎡

Pnear m =

A0 ⎢ Bnear,m ⎣ 0 near,m C0

A1 near ,m B1 near ,m C1

A2 near ,m B2 near ,m C2

⎤

⎥ ⎦.

(44)

Here, m < 40 and the maximum number of nonzeros per row from AFEM is less than 50 on a good quality unstructured grid. Therefore, the preconditioner Pnear m is a highly sparse matrix giving a good approximation to the original hybrid system. The state-of-art parallel direct solver MUMPS [38] is used to compute the inverse of the preconditioner (Pnear or Pnear m ) that is involved in the iterative solution of the preconditioned version of Eq. (28).

712

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

Fig. 2. Illustration of the 3D-1 model which consists of an anomalous prism embedded in a homogeneous background Earth. The frequency of the uniform inducing field is 0.1 Hz. The resistivity of the anomalous prism is ρ = 104 m. The resistivity of the background Earth is 100 m. The profile is arranged along the x-axis at y = 0 m on the flat air–Earth interface.

¯ ¯Ju and ¯J v are available. At a given point in Ω , the electric field After solving the above system of linear equations, E, is computed by 3D interpolation of the terms in Eqs. (25); and the magnetic field is computed through Faraday’s law in Eq. (1). 2.5. Impedance tensor and vertical magnetic transfer function To compute MT and RMT transfer functions from E and H fields, we need to solve the above hybrid surface and volume integral equations (20) and (21) with known surface virtual electric current densities JΓ1 and JΓ0 of at least two different source polarizations. At each source polarization, these known surface virtual electric current densities JΓ1 and JΓ0 are computed from the analytical solutions of magnetic fields of layered Earth models [28] on surfaces Γ1 and Γ0 , respectively. p p These two source polarizations are denoted by symbols E x and H x . p p For the E x source polarization, the primary electric field E is along the x-direction and its amplitude at the air side p of the air–Earth interface on Γ0 is set to 1. For the H x source polarization, the primary magnetic field H p is along the x-direction and its amplitude at the air side of the air–Earth interface on Γ0 is set to 1. Assuming that there is a local coordinate system (u , v , n) at each measuring site (see Fig. 1) on the air–Earth interface, the components of the electric fields and magnetic fields are related by an impedance tensor Z and the vertical magnetic transfer function (VMTF) T [1]. The impedance tensor is determined as follows:

Z uu Z vu

Z uv Zvv



⎡ =⎣

E

p

Eux p

E Evx

H

p

Eu x p

H Ev x

⎤⎡ ⎦⎣

E

p

Hux p

E Hvx

H

p

Hu x p

H Hv x

⎤−1 ⎦

.

(45)

The VMTF is computed from the expression

Tu Tv



⎡ =⎣

E

p

Hux H

E

p

Hu x

p

Hvx H

p

Hv x

⎤−1 ⎡ ⎦

⎣

E

p

Hn x H

p

⎤ ⎦.

(46)

Hn x

Once the impedance tensor is available, the apparent resistivities and phases are extracted from the basic definitions

ρai j =

1

ωμ0

| Z i j |2 ,

for i , j = u , v ,



 φ i j = arctan Im( Z i j )/ Re( Z i j ) .

(47) (48)

3. Numerical experiments 3.1. 3D-1 model The first numerical example is the 3D-1 model [39]. Its geometry is illustrated in Fig. 2. This model has an anomalous prism of resistivity 104 m embedded within a homogeneous background Earth of resistivity 100 m. Both the dielectric permittivity and magnetic permeability of the prism and the Earth are the same as the values in the air space. The incident wave frequency is 0.1 Hz. The uniform inducing fields are vertically incident on the flat air–Earth interface. The mesh has 132 716 edges in Ω1 and 2700 triangles on Γ . Because the performance of the AMFM BEM–FEM algorithm is linearly proportional to the number of iterations of the BICGSTAB solver, we compared the performances of several different preconditioners, which are the block Jacobi preconditioner, the ILU preconditioner, the sparse preconditioner Pnear in Eq. (43) and the highly sparse preconditioner Pnear m in

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

713

Fig. 3. Comparison of performances of different preconditioners used by the BICGSTAB solver to solve the system of linear equations AXt = F (t is the number of iterations) arising from the AMFM BEM–FEM approach on the 3D-1 model in Fig. 2 at a frequency of 0.1 Hz. Due to the similarity of the ILU(0) and block Jacobi preconditioners, they show almost exactly the same convergence behavior. The convergence rates of the Pnear m (m = 20, Eq. (44)) and Pnear (Eq. (43)) preconditioners are superior to those of the other preconditioners.

Table 1 Mesh and execution parameters for the 3D-1 model in Fig. 2 using the standard BEM–FEM, AMFM BEM–FEM and FEM approaches. The domain Ω1 is the inhomogeneous Earth. The air–Earth interface is denoted by Γ . The stopping tolerance for the preconditioned BICGSTAB solver is AXt − F/F < 10−10 (t is the number of iterations). Both the standard BEM–FEM approach and the FEM algorithm employ the direct solver MUMPS [38]. Methods

Edges of Ω1

Triangles of

Γ BEM–FEM AMFM BEM–FEM FEM

132 716 132 716 169 756

2700 2700

Solving AX = F (seconds)

Memory (GB)

Iterations (t)

4601 4227 167

11.0 2.7 2.0

– 6 –

Table 2 near Memory and time cost of solving the preconditioned systems with preconditioners Pnear in the AMFM BEM–FEM approach (see Table 1). The m and P

near model is the 3D-1 model (Fig. 2) at a frequency of 0.1 Hz. The direct solver MUMPS [38] is used to compute the inverses of Pnear in the m and P preconditioned BICGSTAB solver. Due to high variations of the number of nonzeros per row in Pnear , preconditioner Pnear given in Eq. (43) has significantly near poorer performance than preconditioner Pm given in Eq. (44).

Preconditioner

Size

Maximum number of nonzeros per row

Sparsity

Memory (GB)

Time (LU) (seconds)

Total time (BICGSTAB) (seconds)

Pnear m Pnear

138 116 138 116

145 3465

1.50% 2.77%

2 .7 15.7

153 2749

4427 7023

Eq. (44). The block Jacobi preconditioner is constructed by taking the block diagonals of the system matrix A. The ILU preconditioner (here, ILU(0)) is an incomplete LU factorization with 0 level of fill-in of the system matrix A [34]. To construct the highly sparse preconditioner Pnear m , we set m = 20 in Eq. (44) in the experiment. Comparisons of the performances of these preconditioners are shown in Fig. 3. The poor convergence rate of the non-preconditioned BICGSTAB solver implies a large condition number of the original system matrix A. This is because the system matrix is coupled from two physically different discretization systems, where values of the matrix entries arising from the surface integral formula on the air–Earth interface are much smaller in amplitude than those arising from the weak-form finite-element formula in the heterogeneous Earth. The block-Jacobi and ILU preconditioned BICGSTAB solvers have similarly poor convergence behavior owing to the fact that block diagonal structures are only coarse approximations of the partially dense system matrix. In comparison to the un-preconditioned iterative solver, the improvement of the convergence rate is at best marginal. The two sparse preconditioners Pnear and −14 after Pnear m lead both to remarkable convergence rates of the BICGSTAB solver. The relative residual norm is less than 10 near less than 10 iterations. As shown in Table 2, due to the denser structure of preconditioner P with high variations of nonzero entries per row, the memory and time costs for the computation of the inverse of the preconditioner Pnear are

714

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

Fig. 4. Comparison of the apparent resistivities ρa and phases φ of the off-diagonal terms of the impedance tensor and the x-component of the VMTF T x [1] for the 3D-1 model (Fig. 2) at a frequency of 0.1 Hz. The reference FEM solutions are obtained by the goal-oriented edge-based finite-element method [14].

near much higher than those for the highly sparse preconditioner Pnear m . Therefore, the preconditioner Pm is more suitable for our nonuniform unstructured meshes and will be used in the following experiments. The detailed mesh and execution parameters for different methods are listed in Table 1. For this moderate size model, the AMFM BEM–FEM approach outperforms the standard BEM–FEM approach in terms of computational cost. We also compare the speed of the hybrid scheme against the FEM approach in Table 1. For this low frequency problem, the wavelengths of the electromagnetic fields are rather long so a smaller number of unknowns is needed in the FEM approach. Compared to the linear time complexity O ( N e + N e0 ) of the FEM approach (N e is the number of edges in the Earth and N e0 is the number of edges in the air space), the optimal time complexity O ( N e + N tΓ log N tΓ ) (N tΓ is the number of triangle over the air–Earth interface) of the AMFM accelerated BEM–FEM is still less efficient, if N e0 < N tΓ log N tΓ . The numerical solutions of apparent resistivities, phases and VMTFs are shown in Fig. 4. Due to the unavailability of analytical solutions for this model, the numerical solutions are compared to a reference solution computed by a goal-oriented adaptive FEM algorithm [14]. In Fig. 4, the solution of the standard BEM–FEM approach has excellent agreement with those of the AMFM BEM–FEM approach and both of them also have excellent agreement with the reference adaptive FEM solution (in which the air space is included with volume discretization). These consistencies not only verify the idea of the presented hybrid scheme, but also confirm the implementation of the much involved AMFM BEM–FEM algorithm.

3.2. A three-layer model with topography The second example is taken from a simplified 2D resistivity inversion result of an RMT survey in Sweden [12] which consists of three layers having irregular topography as shown in Fig. 5. A profile is arranged along the topographical hill in the y-axis direction between y = −100 m and y = 100 m. The frequency is 200 kHz for vertically incident plane waves. The computation domain is [−800 m, 800 m]3 for the FEM calculation. The underground part of the domain [−800 m, 800 m]2 × [−10 m, 800 m] is used for the hybrid BEM–FEM calculation. As boundary conditions, analytical 1D solutions for the 3-layered Earth model at the left and right edges of the model (Fig. 5) are employed.

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

715

Fig. 5. The y–z section view of the 3-layer irregular interface topography model (its x–z section is identical to its y–z section). The magnetic permeability of the free space is used for all three layers.

Fig. 6. Comparison of the apparent resistivities ρa and phases φ of the off-diagonal terms of the impedance tensor and the y-component of the VMTF (T y ) [1] for the 3-layer irregular topography model (Fig. 5) at a frequency of 200 kHz on a profile along the y-axis. The reference FEM solutions are obtained by the goal-oriented edge finite-element method [14].

Comparisons of apparent resistivities, phases and VMTFs are shown in Fig. 6. There is still good agreement between the AMFM BEM–FEM and the FEM approaches. The average error of the AMFM BEM–FEM solution relative to the FEM solution is less than 2.0%. The largest relative errors occur around changes in topographic slope at y = −5 m and y = 5 m (Fig. 6a) and over the top of the hill (Fig. 6b). The relative error is about 15.0% at y ± 5 m, which is caused by the singular behavior of electromagnetic fields at the corners [14]. The relative error is less than 5.0% over the top of the hill. To generate smoother solutions for the AMFMM BEM–FEM approach, further investigation of applying adaptive mesh refinement in the AMFMM BEM–FEM algorithm is needed [40].

716

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

Table 3 p Execution parameters for the three-layer model with topography in Fig. 5 using the adaptive FEM code [13] and AMFM BEM–FEM for the E x polarization. A fifth-order Gauss quadrature rule is used to evaluate the regular integrals in the system matrix A and the right-hand side vector F. The LU direct solver MUMPS [38] is employed to solve the system of linear equations for the goal-oriented adaptive FEM approach. The maximum number of triangles in a leaf cluster is 20. The cluster tree of the AMFM BEM–FEM approach has 12 levels, 8658 clusters and 6525 leaf clusters. The terminating tolerance for preconditioned BICGSTAB solver is Xt − F/F = 1.8 × 10−06 employing preconditioner Pnear m (m = 20) in Eq. (44). Three-layer model with topography

Unknowns

Memory cost (GB)

Solving AX = F (seconds)

adaptive FEM

1 324 998

>32.0

AMFM BEM–FEM

661 529

15.0

Solvers

Iterations (t)

11 504

LU

–

9779

Pnear m + BICGSTAB

4

The computational cost of both approaches is listed in Table 3. Because a direct solver is used in the FEM code, a large amount of memory is needed which easily approaches the peak memory of our platform. The FEM code has to run in the out-of-core mode. Instead, the hybrid BEM–FEM approach not only shows a very fast convergence rate with a relative residual tolerance of 1.8 × 10−6 in 4 iterations, but also requires less memory (see Table 3). 4. Discussion and conclusions In this study, we introduce a novel hybrid approach for electromagnetic induction modeling which completely eliminates the volume discretization of the air space. The edge-based unstructured finite-element method is used to model the complicated heterogeneous Earth with arbitrary distributions of conductivity, magnetic permeability and dielectric permittivity. A surface integral formula which treats the reduced electric vector potentials as its primary variable is introduced on the air–Earth interface. Due to the surface integral approach, the air–Earth interface can be arbitrarily complicated. The point collocation boundary element method is used to solve the surface integral equation. As a benefit of the hybrid BEM–FEM scheme, possible artificial discretization errors from the air space are almost completely removed, which is critical at high frequency. In other methods that utilize volume grid discretization, very careful mesh design in the air is required due to the small wavelengths of electromagnetic fields at high frequency. To make the hybrid BEM–FEM approach applicable to large-scale electromagnetic induction problems in the Earth, the adaptive multi-level fast multipole method is adopted as a matrix-free accelerator to the multiplication of the system matrix with a guessed solution vector in the framework of iterative solvers such as the BICGSTAB solver. As explicit matrix storage is avoided, the multi-level, fast multipole based hybrid BEM–FEM schemes avoids the O ( N e + ( N tΓ )2 ) memory cost (N e is the number of edges in the Earth, N tΓ is the number of triangle over the air–Earth interface), which is required by standard BEM–FEM schemes. Furthermore, in terms of the idea of multi-level fast multipole methods, the O (( N tΓ )2 ) time cost of the multiplication of the system matrix with the guessed solution vector is reduced to a time cost of O ( N tΓ log N tΓ ). These two significant improvements of reducing memory cost and run time cost enable our newly developed hybrid BEM–FEM scheme to simulate large-scale problems, which previously were only considered possible by the volume differential methods such as FEM and FDM with sparse system matrices. To further improve the performance of hybrid AMFM BEM–FEM schemes, a highly sparse preconditioner is introduced to significantly reduce the iteration numbers of the BICGSTAB solver. This highly sparse preconditioner is extracted from the submatrices of the FEM part and the banded submatrices of the near-field contributions of the BEM part. The experiments show that for both high frequency (200 kHz) and low frequency cases (0.1 Hz), a relative residual of 10−13 can be obtained after less than 20 iterations. Due to the highly sparse system matrices arising from the FEM and FDM approaches, and the nearly linear memory and time cost of solving these highly sparse linear systems, these is no doubt that FEM or FDM are the methods of choice to simulate models with flat or smooth surface topography, especially at low frequency. However, if the surface topography is significant and high frequency is considered, we found the hybrid AMFM BEM–FEM method might be a preferable approach. The reasons are threefold: (1) possible numerical errors from improper volume grid discretization of the free (air) space which is required for FEM and FDM at high frequencies can be avoided; (2) because very dense grids are needed in the free space at high frequencies with other methods, the time and memory costs of the hybrid AMFM BEM–FEM are comparable to or smaller than those of the FEM and FDM algorithms; (3) careful mesh design in the air space is avoided, which is critical at high frequency. Acknowledgements This work was performed in the Exploration and Environmental Geophysics (EEG) Group of ETH Zurich under the overall direction of Prof. Dr. Johan Robertsson. The research was financially supported by the China Scholarship Council Foundation (2008637007), ETH Zurich and the National Natural Science Foundation of China (40874072). References [1] M.N. Berdichevsky, V.I. Dmitriev, Models and Methods of Magnetotellurics, Springer, Berlin and Heidelberg, 2008.

Z. Ren et al. / Journal of Computational Physics 258 (2014) 705–717

717

[2] T. Kalscheuer, M. De LosÁngeles García Juanatey, N. Meqbel, L. Pedersen, Non-linear model error and resolution properties from two-dimensional single and joint inversions of direct current resistivity and radiomagnetotelluric data, Geophys. J. Int. 182 (3) (2010) 1174–1188. [3] R. Mackie, J. Smith, T. Madden, 3D electromagnetic modeling using finite-difference equations – the magnetotelluric example, Radio Sci. 29 (1994) 923–935. [4] G. Newman, D. Alumbaugh, Three-dimensional induction logging problems, Part 2: A finite-difference solution, Geophysics 67 (2002) 484–491. [5] Y. Mitsuhata, T. Uchida, 3D magnetotelluric modeling using the T-Omega finite-element method, Geophysics 69 (2004) 108–119. [6] M. Nam, H. Kim, Y. Song, T. Lee, J. Son, J. Suh, 3D magnetotelluric modelling including surface topography, Geophys. Prospect. 55 (2007) 277–287. [7] A. Franke, R.U. Börner, K. Spitzer, 3D finite element simulation of magnetotelluric fields using unstructured grids, in: Proceedings of the 4th International Symposium on Three-Dimensional Electromagnetics, TU Bergakademie, Freiberg, Germany, 2007, pp. 138–141. [8] Y. Li, K. Key, 2D marine controlled-source electromagnetic modeling: Part 1 – An adaptive finite-element algorithm, Geophysics 72 (2007) WA51–WA62. [9] G. Hohmann, Three-dimensional induced polarization and electromagnetic modeling, Geophysics 40 (1975) 309–324. [10] D. Avdeev, A. Kuvshinov, O. Pankratov, G. Newman, Three-dimensional induction logging problems, Part I: An integral equation solution and model comparisons, Geophysics 67 (2002) 413–426. [11] M. Everett, Theoretical developments in electromagnetic induction geophysics with selected applications in the near surface, Surv. Geophys. 33 (2011) 1–35. [12] T. Kalscheuer, L.B. Pedersen, W. Siripunvaraporn, Radiomagnetotelluric two-dimensional forward and inverse modelling accounting for displacement currents, Geophys. J. Int. 175 (2008) 486–514. [13] Z. Ren, T. Kalscheuer, S. Greenhalgh, H. Maurer, Boundary element solutions for broadband 3D geo-electromagnetic problems accelerated by multi-level fast multipole method, Geophys. J. Int. 192 (2012) 473–499. [14] Z. Ren, T. Kalscheuer, S. Greenhalgh, H. Maurer, A goal-oriented adaptive finite-element approach for plane wave 3D electromagnetic modelling, Geophys. J. Int. 194 (2012) 700–718. [15] Z. Ren, F. Bouillault, A. Razek, A. Bossavit, J. Verite, A new hybrid model using electric field formulation for 3-D eddy current problems, IEEE Trans. Magn. 26 (1990) 470–473. [16] D. White, B. Fasenfest, R. Rieben, M. Stowell, A QR accelerated volume-to-surface boundary condition for the finite-element solution of eddy-current problems, IEEE Trans. Magn. 43 (2007) 1920–1933. [17] X. Sheng, J. Jin, J. Song, C. Lu, W. Chew, On the formulation of hybrid finite-element and boundary-integral methods for 3-D scattering, IEEE Trans. Antennas Propag. 46 (1998) 303–311. [18] M. Becken, O. Ritter, P. Bedrosian, U. Weckmann, Correlation between deep fluids, tremor and creep along the central San Andreas fault, Nature 480 (2011) 87–90. [19] L.B. Pedersen, M. Bastani, L. Dynesius, Some characteristics of the electromagnetic field from radio transmitters in Europe, Geophysics 71 (2006) G279–G284. [20] L. Greengard, The Rapid Evaluation of Potential Fields Particle Systems, MIT Press, 1988. [21] W. Hackbusch, Z. Nowak, On the fast matrix multiplication in the boundary element method by panel clustering, Numer. Math. 54 (1989) 463–491. [22] J. Jin, The Finite Element Method in Electromagnetics, Wiley-IEEE Press, New York, 2002. [23] C. Schwarzbach, Stability of finite element solutions to Maxwell’s equations in frequency domain, PhD thesis, 2009. [24] G. Keller, Rock and mineral properties, in: M.N. Nabighian (Ed.), Electromagnetic Methods in Applied Geophysics Theory, vol. 1, Society of Exploration Geophysicists, Tulsa, 1987, pp. 13–51. [25] O. Biro, K. Preis, Finite element analysis of 3-d eddy currents, IEEE Trans. Magn. 26 (2) (1990) 418–423. [26] T. Onuki, Hybrid finite element and boundary element method applied to electromagnetic problems, IEEE Trans. Magn. 26 (1990) 582–587. [27] C. Carpenter, Comparison of alternative formulations of 3-dimensional magnetic-field and eddy-current problems at power frequencies, IEE Proc B 127 (1980) 332. [28] S. Ward, G. Hohmann, Electromagnetic theory for geophysical applications, in: M.N. Nabighian (Ed.), Electromagnetic Methods in Applied Geophysics Theory, vol. 1, Society of Exploration Geophysicists, Tulsa, 1987, pp. 131–311. [29] J. Stratton, Electromagnetic Theory, Wiley–IEEE Press, 2007. [30] J. Nédélec, A new family of mixed finite elements in R3, Numer. Math. 50 (1986) 57–81. [31] M. Duffy, Quadrature over a pyramid or cube of integrands with a singularity at a vertex, SIAM J. Numer. Anal. 19 (1982) 1260–1262. [32] Y. Liu, Fast Multipole Boundary Element Method, Cambridge University Press, 2009. [33] S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith, H. Zhang, PETSc users manual, Tech. rep, Argonne National Laboratory, 2010. [34] R. Barrett, M. Berry, T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd ed., SIAM, Philadelphia, PA, 1994. [35] M. Benzi, Preconditioning techniques for large linear systems: a survey, J. Comput. Phys. 182 (2) (2002) 418–477. [36] J. Lee, J. Zhang, C. Lu, Incomplete LU preconditioning for large scale dense complex linear systems from electromagnetic wave scattering problems, J. Comput. Phys. 185 (2003) 158–175. [37] Y. Wang, J. Lee, J. Zhang, A short survey on preconditioning techniques for large-scale dense complex linear systems in electromagnetics, Int. J. Comput. Math. 84 (8) (2007) 1211–1223. [38] P. Amestoy, I. Duff, J. L’Excellent, J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM J. Matrix Anal. Appl. 23 (2002) 15–41. [39] M. Zhdanov, I. Varentsov, J. Weaver, N. Golubev, V. Krylov, Methods for modelling electromagnetic fields: Results from COMMEMI – the international project on the comparison of modelling methods for electromagnetic induction, J. Appl. Geophys. 37 (1997) 133–271. [40] E. Stephan, M. Maischak, A posteriori error estimates for FEM–BEM couplings of three-dimensional electromagnetic problems, Comput. Methods Appl. Math. 194 (2) (2005) 441–452.