A Secondary Field Based hp-Finite Element Method for the Simulation of Magnetotelluric Measurements

A Secondary Field Based hp-Finite Element Method for the Simulation of Magnetotelluric Measurements J. Alvarez-Aramberria,b,c,∗, D. Pardob,a,d, H. Barucqc,e aBCAM (Basque Center for Applied Mathematics), Mazarredo 14, 48009, Bilbao, Spain bUniversity of the Basque Country (UPV/EHU), Bilbao, Spain cUniversity of Pau (UPPA), Pau, France dIkerbasque, Bilbao, Spain eInria team-project Magique-3D, Pau, France Abstract In some geophysical problems, it is sometimes possible to divide the subsurface resistivity distribution as a one dimensional (1D) contribution plus some two dimensional (2D) inhomogeneities. Assuming this scenario, we split the electromagnetic fields into their primary and secondary components, the former corresponding to the 1D contribution, and the latter to the 2D inhomogeneities. While the primary field is solved via an analytical solution, for the secondary field we employ a multi-goal oriented self-adaptive hp-Finite Element Method (FEM). To truncate the computational domain, we design a Perfectly Matched Layer (PML) that automatically adapts to high-contrast materials that appear in the subsurface and in the air-ground interface. Numerical results illustrate the robustness of the proposed PML and the gains of the secondary field approach, where we obtain results with comparable accuracy than with a full field based formulation but with a much lower computational cost. Keywords: Finite Element Method (FEM), hp-adaptivity, Magnetotelluric Problem, Secondary Field Formulation, Perfectly Matched Layers (PML). ∗Corresponding author Email address: [email protected] (J. Alvarez-Aramberri) March 13, 2015 This is the accepted manuscript of the article that appeared in final form in Journal of Computational Science 11 : 137-144 (2015), which has been published in final form at https://doi.org/10.1016/j.jocs.2015.02.005. © 2015 Elsevier under CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) 1. Introduction The magnetotelluric (MT) method is a passive exploration technique based on electromagnetic (EM) waves [1, 2, 3]. It aims at estimating the resistivity distribution, and therefore at providing an image of the Earth’s subsurface. MT measurements are governed by Maxwell’s equations with a surface source located at the ionosphere. In particular, when the materials and the source depend only upon two spatial variables, two independent and uncoupled modes are derived, the so-called Transverse Electric (TE) and Transverse Magnetic (TM) polarizations. The solution to the equations arisen from these two modes can be numerically solved with a hp-Finite Element Method (FEM) [4, 5, 6, 7]. With those solutions, it is then possible to compute the impedance and/or the apparent resistivity, two suitable physical quantities to perform the inversion. To correctly capture the complexity of the Earth’s subsurface, we employ adaptive grids, which allow to approximate special features of the solution by refining only in specific areas. To build the refined mesh, we employ a goal-oriented adaptive strategy [8], which minimizes the error of a prescribed quantity of interest represented by a linear functional (see [9, 8, 10, 11, 12] for details). The ability of the goal-oriented algorithm to provide accurate solutions in a region of interest in the context of hp-FEM has been described in various works [13, 14, 15]. The hp-FEM provides exponential convergence rates for elliptic problems with a piecewise analytic solution, whereas hor pversions converge only algebraically. This was proved in 1D by Gui and Babuska [16] and in 2D by Babuska and Gui [17] and Schwab [18].Only the hp-FEM is able to combine small elements (needed to capture geometrical details such as thin edges) with high orders of approximation (necessary to decrease the dispersion error for wave propagation problems [19, 20, 21]). Besides, it is robust for singularly perturbed problems, that is, it still performs appropriately when a parameter involved in a given elliptic problem approaches a critical value [18]. In some geophysical applications, as in MTs, the data is acquired at several receivers located at the Earth’s surface. It becomes then necessary to obtain accurate results at multiple positions, being this the reason to extend the goal-oriented strategy to a multigoal-oriented one. There exist two possible approaches towards multigoal-oriented adaptivity. The first one consists of using one grid for each goal, as in [22], where the implementation needs to handle multiple grids, which in general may be complicated. The second 2 one consists of defining a new quantity of interest that takes into account all goals (see [23, 24]). Based on this second approach, we implement the algorithm proposed by Pardo in [25]. In geophysics in general and in MT in particular, when the subsurface distribution of the resistivity depends upon multiple spatial variables, it is sometimes possible to interpret it as a 1D formation plus some 2D (or 3D) heterogeneities. In this work, we consider a horizontally layered Earth model with 2D heterogeneities. Then, in order to solve the TE and TM modes, we split the electric and magnetic fields into their primary and secondary components. The first corresponds to the fields arisen from some reference conductivity model (1D), while the second arises from the difference between the actual conductivity distribution with respect to the reference model (2D). Since the 1D solution is known analytically, the main advantage of this approach is that we only need to accurately solve the secondary field variations (the term “secondary field” is also known as “scattered field” in the electrical engineering community), which in general are easier to solve, since they exhibit less variations (smaller gradients) than the primary field. Hence, it is generally possible to employ coarser grids, and hence reduce the computational cost. Additionally, in MTs the computational domains are usually very large if one models the incident plane wave source. In the secondary field formulation, the source term is not at the ionosphere, but where the inhomogeneities are considered. Since it is not necessary anymore to model the ionosphere source, this allows us to considerably reduce the computational domain. Finally, since we separate the primary from the secondary field, we may obtain additional physical relevant information by analyzing each field (primary and secondary) separately. The main contribution of this work is then to solve, via the hp-FEM, the MT direct problem using the secondary field formulation to simulate MT measurements. The mentioned benefits of this approach will then be notorious in the inversion. On the one hand, there exists the possibility of analyzing 1D and 2D effects separately. On the other hand, since the solution of the inverse problem is based on reiterated solutions of the direct problem, reducing the computational cost of solving the direct problem produces large savings in the computational costs. Additionally, we provide an automatic technique to truncate the computational domain, a further problem that appears when applying a FEM to unbounded region problems such as MT. Different approaches can be em3 ployed for this purpose. We employ a Perfectly Matched Layer (PML), which is an exact method at the continuous level, and thus, it matches the highaccuracy delivered by the hp-adaptive FEM. The work of Gomez-Revuelto et al. [26] shows the suitability of the utilization of PMLs in this context. PMLs were proposed by Berenger [27] (1994) in an electromagnetic context as an artificial layer intended to reduce reflections from the boundary of a truncated computational domain. In this method, one has to select the decay profile of the wave into the PML region. This profile needs to ensure that reflections from the boundary are arbitrarily small, which implies that the solution decays arbitrarily fast, creating then a “boundary layer” with strong gradients within the PML region. Thus, while a low decay produces incoming waves reflected from the boundary, an excessive decay requires a very fine grid to approximate the boundary layer. To find an equilibrium between a fast and a slow decay, it is necessary to properly adjust the PML parameters, which is usually tricky since they depend on the problem itself. Moreover, when we have a layered material with high contrasts on the material properties, this selection of the parameters is even more challenging. Thus, in this work we also provide a method to automatically adjust the PML parameters, even in the most complex scenarios where the material contrast properties among neighboring materials are as high as sixteen orders of magnitude. These type of scenarios often appear in geophysical electromagnetic (EM) applications that involve both, air and ground. We show that the proposed PML produces an appropriate decay of the solution in the air and in the subsurface without introducing spurious reflections, and thus, providing accurate solutions. The present work is organized as follows. In Section 2 we define the formulation of the problem. Section 3 describes the formulation of the PML and how the parameters are adjust in the Automatically Adapted PML. We derive the secondary field formulation in Section 4 and numerical results based on the MT problem are illustrated in Section 5. Section 6 is devoted to the conclusions. 2. Formulation of the Method MT measurements are governed by the electromagnetic phenomena, which is described by Maxwell’s equations. Assuming a time-harmonic dependence 4 of the form ejωt, these equations can be expressed in frequency domain as:    ∇×E=−jωµH −Mimp (Faraday), ∇×H= (σ+jωε)E+Jimp (Amp`ere), (1) where Eand Hare the electric and magnetic fields, respectively. These fields are driven by an impressed prescribed electric and magnetic density current sources, Jimp = (0, Jy,0) and Mimp = (0, My,0), respectively. We emphasize that, as explained in [28], magnetic impressed currents are only mathematical symbols utilized to represent sources. jis the imaginary unit, ωis the angular frequency, and σstands for the conductivity of the media. We assume that σ=   σ0 0 0σ0 0 0 σ  ,ε=   ε0 0 0ε0 0 0 ε  ,µ=   µ0 0 0µ0 0 0 µ  ,(2) where the electrical permitivity εand the magnetic permeability µare assumed to be that of the vacuum (ε0and µ0respectively) and σ(x, y, z) to be piecewise constant, non-negative, and bounded above. 2.1. Transverse Electric (TE) and Transverse Magnetic (TM) Modes Pre-multiplying both sides of Faraday’s Law by µ−1, applying the curl, and using Amp`ere’s Law, we obtain the reduced wave equation, ∇×(µ−1∇×E)−k2E=−jωJimp −∇×µ−1Mimp,(3) where k2=ω2ε−jωσ. A similar equation is obtained in terms of the magnetic field by multiplying both sides of Amp`ere’s Law by ˆ σ−1= (σ+ jωε)−1and applying the curl to Faraday’s Law ∇×(ˆ σ−1∇×H) + jωµH =−Mimp +∇×ˆ σ−1Jimp.(4) When the materials and the source depend only upon two spatial variables (x, z), then ∂/∂y = 0 and two independent and uncoupled modes are derived from Maxwell’s equations. The uncoupled TE mode involves (Ey, Hx, Hz) field components, while TM only considers (Hy, Ex, Ez). Our aim is to find the ycomponent of the electric and magnetic fields Ey(x, z), Hy(x, z)∈ H1(Ω) that satisfy the BCs and equations (3), and (4), respectively. 5 To derive the variational formulation, we first define the L2-inner product of two possible complex and vector valued functions g1and g2as: hg1,g2iL2(Ω) =ZΩ g1 ∗g2dΩ,(5) where g∗denotes the adjoint (transpose of the complex conjugate) of g. 2.1.1. TE Variational Formulation To obtain the corresponding variational formulation, we pre-multiply (3) by the complex conjugate of a test function F∈V(Ω), where V(Ω) = H1 ΓD(Ω) = {F∈L2(Ω) : F|ΓD= 0,∇F∈L2(Ω)}is the space of admissible test functions. Then, we integrate by parts and we incorporate the homogeneous Dirichlet BC (the ones considered in the present work) over ΓD=∂Ω. Thus, we obtain:    Find Ey∈V(Ω),such that: h∇F, µ−1∇EyiL2(Ω) − hF, k2EyiL2(Ω) =−jωhF, Jimp yiL2(Ω) ∀F∈V(Ω), (6) 2.1.2. TM Variational Formulation In a similar way, from (4) we obtain the corresponding variational formulation for the magnetic field    Find Hy∈V(Ω),such that: h∇F, ˆσ−1∇HyiL2(Ω) +jωhF, µHyiL2(Ω) =−hF, Mimp yiL2(Ω) ∀F∈V(Ω), (7) We employ an hp-Finite Element Method [4] to solve both problems (6) and (7). The objective of the traditional goal-oriented method is to construct an optimal hp-grid in the sense that it minimizes the problem size needed to achieve a given tolerance error for a given quantity of interest (solution at the receiver) Li(u), being ueither Eyor Hy. This quantity is a linear and continuous functional [14, 15] in uassociated to the i-th receiver and defined as: Li(u) = 1 |ΩRi|ZΩRi u dΩ,(8) where ΩRiis the domain occupied by the i-th receiver. 6 Since we have more than one receiver, we need to properly compute several quantities of interest. Therefore, we employ a multigoal-oriented strategy, proposed in [25], where a new linear quantity of interest that takes into account all receivers is employed. From the solution to the variational problems, we compute the impedance and/or the apparent resistivity, which are two post-processed transfer functions that are typically used during inversion in MT problems. The impedance Zis defined as Zi T E =Zi yx =Li(Ey) Li(Hx),Zi T M =Zi xy =Li(Ex) Li(Hy),(9) where Hxand Eyare obtained from Maxwell’s equations as Hx=1 jωµ ∂Ey ∂z , Ex=−1 σ+jωε ∂Hy ∂z .(10) The apparent resistivity ρapp is defined as ρapp mn =|Zmn|2 ωµ .(11) For the sake of simplicity in the notation, we omit the ysubscript from Ey, Hy,Jimp y, and Mimp yfrom now on. 3. Truncation of the Domain When applying a FEM to unbounded region problems such as MT, the computational domain must be truncated. We employ PMLs for this purpose, and we follow the interpretation introduced by Teixera and Chew in [29, 30], where they consider a PML as an analytic continuation of the governing equations into the complex plane (see also [31]). PMLs transform propagating and evanescent waves into exponentially fast decaying evanescent waves. Since waves are strongly attenuated inside the PML region, the bounded computational domain can be limited by a surface on which one can set perfectly reflecting BCs (in our case, we set homogeneous Dirichlet BCs). Indeed, any reflected wave is so much absorbed inside the layer that it does not pollute the solution inside the domain of study. Then, the selected BCs for both problems imply that the tangential component of the fields are zero on the outer part of the boundary. For a recent review of the state of the art of this truncation technique, see [32] and [33]. 7 3.1. PML Definition Let the Cartesian coordinate system x= (x, z) be the reference system of coordinates in a 2D scenario, where for simplicity, we select the vertical coordinate zas the direction perpendicular to the ground-air interface. Given an arbitrary complex system of coordinates ζ= (ζ1, ζ2), we define our change of coordinates by x=ψ(ζ)and we denote the Jacobian matrix and its determinant by Jand det(J). The change of coordinates is assumed to be represented by an injective differentiable function with continuous partial derivatives and nonzero determinant at any point. We define a one dimensional change of variables in the positive direction of the i-th coordinate as ζxi(xi) = Zxi 0h(η)dη, for i= 1,2, x1=x, x2=z, (12) where h(·) is a possibly complex valued function to be determined in section 3.3. The case corresponding to the negative direction can be defined analogously. The Jacobian is given by [J]i,j ="∂ζi ∂xj#i,j , for i, j = 1,2. Thus, it is expressed as J=   h(x) 0 0h(z)   ,where det(J) = h(x)h(z),(13) denotes the determinant of the Jacobian. With this particular change of coordinates, the Jacobian is diagonal. However, we derive herein the variational formulation for a general, non orthogonal change of variables. This is useful for other purposes, e.g. development of non-orthogonal Fourier FEMs in certain geometries (see [34, 35]). 3.2. Variational Formulation in an Arbitrary System of Coordinates We define the change of coordinates e E:= E◦ψ=e E(ζ), e F:= F◦ ψ=e F(ζ), and e Jimp := Jimp ◦ψ=Jimp(ζ). Using Einstein’s summation convention, according to the chain rule, denoting with the upper bar the complex conjugate, and taking into account that if f∈C1(Ω), then for all i ∂f ∂ζi =∂f ∂ζi ,∂f ∂ζi =∂f ∂ζi ,(14) 8 we obtain that ∇ζe E=∂e E ∂xi ∂xi ∂ζn exn= (J−1)T∇E, ∇ζe E=∂e E ∂xi ∂xi ∂ζn exn= (J−1)∗∇E. (15) Therefore, multiplying (3) by the complex conjugate of a test function e F, integrating by parts, and incorporating the homogeneous Dirichlet BC over e ΓD, we obtain: h∇ζe F, e µ−1∇ζe EiL2(e Ω) =h(J−1)∗∇F, µ−1(J−1)T∇E det(J)iL2(Ω) = =h∇F, J−1µ−1(J−1)T∇E det(J)iL2(Ω) , he F, e k2e EiL2(e Ω) =hF, k2E det(J)iL2(Ω), he F, e JimpiL2(e Ω) =hF, Jimpdet(J)iL2(Ω), (16) where e µ:= µ◦ψ,e k:= k◦ψ,e Ω := Ω ◦ψ, and e ΓD:= ΓD◦ψ. Following the ideas of [36] concerning the inclusion of metric-dependent variables within material coefficients, we define the following functions: µT E NEW =JTµJ1 det(J)=       µh(x) h(z)0 0µh(z) h(x)        , k2 NEW =k2det(J) = k2h(x)h(z), Jimp NEW =Jimp det(J) = Jimph(x)h(z). (17) The new source and new material tensors incorporate the information about the change of coordinates. Thus, the variational formulation can be expressed in terms of an arbitrary system of coordinates by simply considering the new source and materials. The new variational formulation for the electric field 9 σ1σ2σ3 Model 1 1 1 1 Model 2 1 1/10 1/3 Model 3 1 1/10 1/10 Model 4 1 1/100 1/3 Table 1: Different models for the formation of the subsurface. Conductivities are given in S/m. frequencies between the numerical hp-FEM solution and the exact solution. We obtain relative errors below 1.5%, a superb accuracy for these type of simulations. 10−5 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 Frequency (Hz.) Relative error in percent Model 1 Model 2 Model 3 Model 4 Figure 5: Relative error between the exact and numerical solutions for different subsurface formations against frequency. To study the behavior of the solution into the PML region, we consider Model 4 with frequency equal to 10−4Hz and we display the logarithm of the module of the impedance along all sides of the computational domain. Thus, we represent log(|ZT E |) in Figure 6. We appreciate that the PML behaves properly everywhere, with a smooth decay for the solution and without introducing numerical reflections even in the areas with high contrast between material properties. Panels (a) and (b) correspond to the intersection between air and ground. There, the contrast between resistivities is about sixteen orders of magnitude and even in this scenario, the decay seems to be superb. 16 (a) (b) (c) (d) Figure 6: log(|ZT E |) corresponding to Model 4, with a 5 km thick PML and α= 10−5. Panel (a) corresponds to the left side of the domain. (b), (c), and (d) correspond to the right, top and bottom parts of the domain, respectively. The black line indicates the region where the PML starts. 5.3. Secondary Field Formulation We consider now a 2D scenario with the following conductivity distribution: σ1= 1/3, σ2= 1/2, σ3= 1/4, σ4= 1/200 S/m. The most sensitive frequency to the target area, that is, the frequency at which the presence of the target affects most the measurements at the receivers, corresponds to 0.05 Hz. Figure 7 shows the final grids after executing the multi-goal oriented adaptivity for the full formulation (left) and for the secondary field based problem (right) at this frequency. The left panel shows a zoom of the final grid with the origin of coordinates at the center. The size of the represented domain is of 40 ×70 km2. The grid in the right panel is the complete grid for the secondary field problem (50 ×70) km2. Figure 8 displays the relative errors in the apparent resistivity between the full field and secondary field solutions. There, positions 1 to 4 correspond to measurements obtained at 0, 4, 8 and 20 km from the center of the domain, respectively. Due to the low errors observed in Figure 8, we conclude that both approaches provide analogous results. However, the number of unknowns needed to achieve these small errors are not the same. We now consider the same model with the same frequency of 0.05 Hz. We compute an overkill solution with a much finer grid obtained after performing 17 Figure 7: Final multi-goal oriented hp-grids. Different colors indicate different values of p. Left: full formulation based problem (zoom). Right: secondary field based problem. adaptivity. We use it to estimate the relative errors corresponding to the secondary field and full formulations after several hand/or pglobal refinements. Figure 9 displays the results of these computations for the receiver located at the center of the domain. We appreciate that, for instance, to achive a (small) relative error of 0.1%, one only requires around 7000 unknowns with the secondary field formulation, while to solve the full formulation problem with the same accuracy, we need around 17000 unknowns. Therefore, with the first approach we only need approximately 40% of the unknowns. 6. Conclusions The multi-goal oriented hp-FEM provides accurate solutions for the MT problem at different receivers simultaneously. We show that by employing the secondary field approach, we obtain significant benefits in comparison with directly using the full field formulation: we can obtain additional physical relevant information by analyzing each field (primary and secondary) separately, and furthermore, this is obtained employing a significantly lower number of unknowns. Since the solution of the inverse problem is based on iterated solutions of the direct problem, reducing the computational cost of solving the direct problem induces high savings in the inversion process. We also provide a method to automatically truncate the computational domain employing PMLs. To find an equilibrium between a fast and a slow decay on the PML region is usually tricky. It depends on the problem itself and it is even more complicated when there exists high contrasts on adjacent 18 10−5 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 Frequency (Hz.) Relative error in percent Position 1 Position 2 Position 3 Position 4 Figure 8: Comparison between the results when using the full formulation and the secondary field formulation. material properties. We have shown that in these complicated scenarios, the Automatically Adapted PML provides an adequate decay, not so fast to require a too fine grid and not so slow to introduce artificial reflections. Since the choice of PML parameters is automatic, the proposed approach is also suitable for inverse problems. Even if the reduction of the computational cost is itself beneficial, the main advantage of solving the inverse problem with this approach consists on the fact that it allows for separate analysis of 1D and 2D effects. This will be analyzed in future research. Acknowledgments Julen Alvarez-Aramberri and David Pardo were partially funded by the Project of the Spanish Ministry of Economy and Competitiveness with reference MTM2013-40824-P, the BCAM “Severo Ochoa” accreditation of excellence SEV-2013-0323, the CYTED 2011 project 712RT0449, and the Basque Government through the BERC 2014-2017 program and the Consolidated Research Group Grant IT649-13 on “Mathematical Modeling, Simulation, and Industrial Applications (M2SI)”. David Pardo has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 644602, by the RISE Horizon 2020 European Project GEAGAM (644602). Julen Alvarez-Aramberri was also partially funded by the University of the Basque Country UPV/EHU under the grant PIFG05/2011. 19 0 0.5 1 1.5 2 x 104 0,01 0,1 1 10 Number of degrees of freedom Relative error in percent secondary full Figure 9: Relative error in logarithmic scale for the apparent resistivity computed with the full and secondary field formulations. References [1] L. Cagniard, Basic theory of the magneto-telluric method of geophysical prospecting, Geophysics 18 (3) (1953) 605–635. [2] K. Vozoff, The magnetotelluric method in the exploration of sedimentary basins, Geophysics 37 (1) (1972) 98–141. [3] F. Simpson, K. Bahr, Practical magnetotellurics, Cambridge University Press, 2005. [4] L. Demkowicz, Computing with hp-Adaptive Finite Elements: Volume 1. One and Two Dimensional Elliptic and Maxwell problems, CRC Press, 2006. [5] I. Gomez-Revuelto, L. Garcia-Castillo, S. Llorente-Romano, D. Pardo, 3d hp-adaptive finite element simulations of bend, step, and magic-T electromagnetic waveguide structures, Journal of Computational Science 5 (2) (2014) 65 – 75. [6] J. Alvarez-Aramberri, D. Pardo, H. Barucq, Inversion of magnetotelluric measurements using multigoal oriented hp-adaptivity, Procedia Computer Science 18 (2013) 1564–1573. 20 [7] A. Szymczak, A. Paszy´nska, M. Paszy´nski, D. Pardo, Preventing deadlock during anisotropic 2d mesh adaptation in hp-adaptive FEM, Journal of Computational Science 4 (3) (2013) 170–179. [8] S. Prudhomme, J. Oden, On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors, Computer Methods in Applied Mechanics and Engineering 176 (1) (1999) 313–331. [9] J. T. Oden, S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method, Computers & Mathematics with Applications 41 (5) (2001) 735–756. [10] M. Paraschivoiu, A. T. Patera, A hierarchical duality approach to bounds for the outputs of partial differential equations, Computer Methods in Applied Mechanics and Engineering 158 (3) (1998) 389–407. [11] R. Rannacher, F.-T. Suttmeier, A posteriori error control in finite element methods via duality techniques: Application to perfect plasticity, Computational Mechanics 21 (2) (1998) 123–133. [12] V. Heuveline, R. Rannacher, Duality-based adaptivity in the hp-finite element method, Journal of Numerical Mathematics JNMA 11 (2) (2003) 95–113. [13] P. ˇ Solın, L. Demkowicz, Goal-oriented hp-adaptivity for elliptic problems, Computer Methods in Applied Mechanics and Engineering 193 (6) (2004) 449–468. [14] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, L. Tabarovsky, A goaloriented hp-adaptive finite element method with electromagnetic applications. Part i: electrostatics, International Journal for Numerical Methods in Engineering 65 (8) (2006) 1269–1309. [15] D. Pardo, L. Demkowicz, C. Torres-Verdin, M. Paszynski, A selfadaptive goal-oriented hp-finite element method with electromagnetic applications. Part ii: electrodynamics, Computer methods in Applied Mechanics and Engineering 196 (37) (2007) 3585–3597. [16] W. Gui, I. Babuˇska, The h, p and hp versions of the finite element method in 1 dimension. i-iii, Numerische Mathematik 49 (6) (1986) 577– 683. 21 [17] I. Babuˇska, B. Guo, Approximation properties of the h-p version of the finite element method, Computer Methods in Applied Mechanics and Engineering 133 (3) (1996) 319–346. [18] C. Schwab, p-and hp-finite element methods: Theory and Applications in Solid and Fluid Mechanics, Clarendon Press Oxford, 1998. [19] F. Ihlenburg, I. Babuˇska, Finite element solution of the Helmholtz equation with high wave number part i: The h-version of the FEM, Computers & Mathematics with Applications 30 (9) (1995) 9–37. [20] F. Ihlenburg, I. Babuska, Finite element solution of the Helmholtz equation with high wave number part ii: The hp version of the FEM, SIAM Journal on Numerical Analysis 34 (1) (1997) 315–358. [21] J. Oden, F. Yusheng, Local and pollution error estimation for finite element approximations of elliptic boundary value problems, Journal of Computational and Applied Mathematics 74 (1) (1996) 245–293. [22] P. Solin, L. Dubcova, J. Cerveny, I. Dolezel, Adaptive hp-FEM with arbitrary-level hanging nodes for Maxwell’s equations, Adv. Appl. Math. Mech 2 (4) (2010) 518–532. [23] R. Hartmann, P. Houston, Goal-oriented a posteriori error estimation for multiple target functionals, Springer, 2003. [24] R. Hartmann, Multitarget error estimation and adaptivity in aerodynamic flow simulations, SIAM Journal on Scientific Computing 31 (1) (2008) 708–731. [25] D. Pardo, Multigoal-oriented adaptivity for hp-finite element methods, Procedia Computer Science 1 (1) (2010) 1953–1961. [26] I. Gomez-Revuelto, L. E. Garcia-Castillo, L. F. Demkowicz, A comparison between pml, infinite elements and an iterative bem as mesh truncation methods for hp self-adaptive procedures in electromagnetics, Progress In Electromagnetics Research 126 (2012) 499–519. [27] J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics 114 (2) (1994) 185–200. 22 [28] R. Harrington, J. Harrington, Field computation by moment methods, Oxford University Press, 1996. [29] W. Chew, W. Weedon, A 3d perfectly matched medium from modified Maxwell’s equations with stretched coordinates, Microwave and Optical Technology Letters 7 (13) (1994) 599–604. [30] F. Teixeira, W. Chew, Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves, Microwave and Optical Technology Letters 17 (4) (1998) 231–236. [31] F. Teixeira, W. Chew, Pml-fdtd in cylindrical and spherical grids, Microwave and Guided Wave Letters, IEEE 7 (9) (1997) 285–287. [32] A. Berm´udez, L. Hervella-Nieto, A. Prieto, R. Rodr´ıguez, Perfectly matched layers for time-harmonic second order elliptic problems, Archives of Computational Methods in Engineering 17 (1) (2010) 77– 107. [33] P. Joly, An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation, SeMA Journal 57 (1) (2012) 5–48. [34] D. Pardo, C. Torres-Verd´ın, M. Nam, M. Paszynski, V. Calo, Fourier series expansion in a non-orthogonal system of coordinates for the simulation of 3d alternating current borehole resistivity measurements, Computer Methods in Applied Mechanics and Engineering 197 (45) (2008) 3836–3849. [35] D. Pardo, V. Calo, C. Torres-Verdin, M. Nam, Fourier series expansion in a non-orthogonal system of coordinates for the simulation of 3d-dc borehole resistivity measurements, Computer Methods in Applied Mechanics and Engineering 197 (21) (2008) 1906–1925. [36] A. Ward, J. Pendry, Calculating photonic green’s functions using a nonorthogonal finite-difference time-domain method, Physical Review B 58 (11) (1998) 7252. [37] S. G. Johnson, Notes on perfectly matched layers (pmls), Lecture notes, Massachusetts Institute of Technology, Massachusetts. 23 [38] W. C. Chew, Waves and fields in inhomogeneous media, Vol. 522, IEEE press New York, 1995. 24