A Self-Adaptive Goal-Oriented hp Finite Element Method with Electromagnetic Applications. Part II: Electrodynamics.

D. Pardo a,b, L. Demkowicz a, C. Torres-Verd´ın a,b, M. Paszynski a,c a Institute

for Computational Engineering and Sciences (ICES), The University of Texas at Austin, Austin TX 78712

b Department c On

of Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin TX 78712

leave from AGH University of Science and Technology, Department of Computer Methods in Metallurgy, Cracow, Poland

Abstract We present the formulation, implementation, and applications of a self-adaptive, goal-oriented, hp-Finite Element (FE) Method for Electromagnetic (EM) problems. The algorithm delivers (without any user interaction) a sequence of optimal hpgrids. This sequence of grids minimizes the error in a prescribed quantity of interest with respect to the problem size, and it converges exponentially in terms of the relative error in a user-prescribed quantity of interest against the CPU time, including problems involving high material contrasts, boundary layers, and/or several singularities. The goal-oriented refinement strategy is an extension of a fully automatic, energy-norm based, hp-adaptive algorithm. We illustrate the efficiency of the method with 2D numerical simulations of Maxwell’s equations using both continuous and edge elements. Applications include alternate current (AC) resistivity logging instruments in a borehole environment with steel casing for the assessment of rock formation properties behind casing. Logging instruments, steel casing, and rock formation properties are assumed to exhibit axial symmetry around the axis of a vertical borehole. For the presented challenging class of problems, the self-adaptive goal-oriented hp-FEM delivers results with 5-7 digits of accuracy in the quantity-of-interest. Key words: hp-finite elements, exponential convergence, goal-oriented adaptivity, computational electromagnetism, Maxwell’s equations, through casing resistivity tools (TCRT).

Preprint submitted to Elsevier Science

18 November 2005

1

Introduction

During the last decades, different algorithms have been designed and implemented to generate optimal grids in the solution of relevant engineering problems. Among those algorithms, a self-adaptive, energy-norm based, hp-Finite Element (FE) refinement strategy has been developed at the Institute for Computational Engineering and Sciences (ICES) of The University of Texas at Austin. The strategy produces automatically a sequence of hp-meshes that delivers exponential convergence rates in terms of the energy-norm error against the number of unknowns (as well as the CPU time), independently of the number and type of singularities in the problem. Thus, it provides high accuracy approximations of solutions corresponding to a variety of engineering applications. Furthermore, the self-adaptive strategy is problem independent, and it can be applied to FE discretizations of H 1 -, H(curl)-, and H(div)-spaces, as well as to nonlinear problems (see [7,17] for details). The energy-norm is a quantity of limited relevance for most engineering applications, especially when a particular objective is pursued, for instance, to simulate the electromagnetic (EM) response of geophysical resistivity logging instruments in a borehole environment. In these instruments, the amplitude of the measurement (for example, the electric field) is typically several orders of magnitude smaller at the receiver antennas than at the transmitter antennas. Thus, small relative errors of the solution in the energy-norm do not imply small relative errors of the solution at the receiver antennas. Indeed, it is not uncommon to construct energy-norm based adaptive grids delivering a relative error in the energy-norm below 1%, while the solution at the receiver antennas still exhibits a relative error above 1000% (see [13]). Consequently, a self-adaptive strategy is needed to approximate a specific feature of the solution. Refinement strategies of this type are called goaloriented adaptive algorithms [11,16], and are based on minimizing the error of a prescribed quantity of interest mathematically expressed in terms of a linear functional (see [2,9,12,11,16,18] for details). This paper formulates, implements, and studies (both theoretically and numerically) a self-adaptive hp goal-oriented algorithm intended to solve electrodynamic problems. The algorithm is an extension of the fully automatic (energy-norm based) hp-adaptive strategy described in [7,17], and a continuation of concepts presented in [14,19] for elliptic problems. The organization of this document is as follows: in Section 2, we introduce Maxwell’s equations, governing the EM phenomena. We also derive four corresponding variational formulations for axisymmetric problems. A self-adaptive goal-oriented hp-algorithm for electrodynamic problems is described in Sec2

tion 3. The corresponding details of implementation are discussed in the same section. In Section 4, we illustrate the efficiency of the method with axisymmetric 2D numerical simulations of TCRT. Finally, in Section 5 we draw the main conclusions, and outline future lines of research.

2

Maxwell’s Equations

In this section, we first introduce the time-harmonic Maxwell’s equations in the frequency domain. They form a set of first-order Partial Differential Equations (PDE’s). Then, we describe boundary conditions needed for the simulation of our logging applications of interest. Finally, we derive a variational formulation in terms of either the electric or the magnetic field, and we reduce the dimension of the computational problem by considering axial symmetry.

2.1 Time-harmonic Maxwell’s equations √ Assuming a time-harmonic dependence of the form ejωt , where j = −1 is the imaginary unit, t denotes time, and ω 6= 0 is the angular frequency, Maxwell’s equations can be written as    ∇×H            ∇×E

    ∇ · (¯ ǫE)          ∇ · (µH) ¯

¯ + jω¯ = (σ ǫ)E + Jimp

Ampere’s Law,

¯ H − Mimp = −jω µ

Faraday’s Law,

=ρ

Gauss’ Law of Electricity, and

=0

Gauss’ Law of Magnetism.

(2.1)

Here H and E denote the magnetic and electric field, respectively, Jimp is a prescribed, impressed electric current density, Mimp is a prescribed, impressed ¯ and σ ¯ stand for dielectric permittivity, magnetic current density, tensors ¯ ǫ, µ, magnetic permeability, and electrical conductivity of the medium, respectively, ¯ 6= 0, and and ρ denotes the electric charge distribution. We assume det(µ) ¯ ¯ det(σ + jω ǫ) 6= 0. The equations described in (2.1) are to be understood in the distributional sense, i.e. they are satisfied in the classical sense in subdomains of regular material data, and they also imply appropriate interface conditions across material interfaces. 3

Energy considerations lead to the assumption that both electric field E and magnetic field H must be square integrable. According to equations (2.1)2 and (2.1)4 , Mimp must be divergence free. Maxwell’s equations are not independent. Taking the divergence of Faraday’s Law yields Gauss’ Law of magnetism. By taking the divergence of Ampere’s Law, and by utilizing Gauss’ Electric Law we arrive at the so called continuity equation, ¯ E) + jωρ + ∇ · Jimp = 0 . ∇ · (σ

(2.2)

2.2 Boundary Conditions (BC’s) There exists a variety of BC’s that can be incorporated into Maxwell’s equations. In the following, we describe those BC’s that are of interest for the logging applications discussed in this paper. At this point, we are considering general 3D domains. A discussion on boundary terms corresponding to the axisymmetry condition is included in Section 2.4.

2.2.1 Perfect Electric Conductor (PEC) Maxwell’s equations are to be satisfied in the whole space minus domains occupied by a PEC. A PEC is an idealization of a highly conductive media. ¯ = diag(σc , σc , σc ) with σc → ∞, the corresponding Inside a region where σ electric field converges to zero 1 by applying Ampere’s law. Faraday’s law implies that the tangential component of the electric field E must remain continuous across material interfaces in the absence of impressed magnetic surface currents. Consequently, the tangential component of the electric field must vanish along the PEC boundary, i.e., n ˆ×E = 0 ,

(2.3)

where n ˆ is the unit normal (outward) vector. Since the electric field vanishes inside a PEC, Faraday’s law implies that the magnetic field should also vanish inside a PEC in the absence of magnetic currents. The same Faraday’s law implies that the normal component of the magnetic field premultiplied by the permeability tensor must remain continu1

This result is true under the physical consideration that impressed volume cur¯ E should remain finite, i.e., hJimp , ψi, hσ ¯ E, ψi < ∞ for every test rent Jimp and σ function ψ. See [15] for details.

4

¯ must ous across material interfaces. Therefore, the normal component of µH vanish along the PEC boundary, i.e., ¯ n ˆ · (µH) =0.

(2.4)

The tangential component of the magnetic field (surface current) and normal component of the electric field (surface charge density) need not be zero, and may be determined a-posteriori.

2.2.2 Source Antennas Antennas are modeled by prescribing an impressed volume electric current Jimp or an impressed volume magnetic current Mimp . Using the equivalence principle (see, for example, [8]), we can replace the original impressed electric volume current Jimp with an equivalent electric surface current Jimp = [ˆ n×H]S , S

(2.5)

defined on an arbitrary surface S enclosing the support of Jimp , where [ˆ n×H]S denotes • the jump of n ˆ×H across S in the case of an interface condition, or • simply n ˆ×H on S in the case of a boundary condition. Similarly, an impressed magnetic volume current Mimp can be replaced by the equivalent magnetic surface current Mimp = −[ˆ n×E]S , S

(2.6)

defined on an arbitrary surface S enclosing the support of Mimp .

2.2.3 Closure of the Domain We consider a bounded computational domain Ω. A variety of BC’s can be imposed on the boundary ∂Ω such that the difference between the solution of such a problem and the solution of the original problem defined over R 3 is small. For example, it is possible to use an infinite element technique (as described in [3]). Also, since the EM fields and their derivatives decay exponentially in the presence of lossy media (non-zero conductivity), we may simply impose a homogeneous Dirichlet or Neumann BC on the boundary of a sufficiently large computational domain. 5

In the field of geophysical logging applications, it is customary to impose a homogeneous Dirichlet BC on the boundary of a large computational domain (for example, 2-20 meters in each direction from a 2 Mhz source antenna in the presence of a resistive media). We will follow the same approach. According to the boundary conditions discussed above, we will represent the boundary Γ = ∂Ω as the disjoint union of • ΓE , where Mimp n×E]ΓE (with Mimp ΓE = −[ˆ ΓE possibly zero), and

(2.7)

• ΓH , where Jimp n×H]ΓH (with Jimp ΓH = [ˆ ΓH possibly zero).

(2.8)

2.3 Variational Formulation From Maxwell’s equations and the BC’s described above, we derive the corresponding standard variational formulation in terms of the electric or magnetic field as follows. First, we notice from Faraday’s law that ∇×E ∈ (L2 (Ω))3 if and only if Mimp ∈ (L2 (Ω))3 . Since our objective is to find a solution E ∈ H(curl; Ω) = {F ∈ (L2 (Ω))3 : ∇×F ∈ (L2 (Ω))3 }, we shall assume in the case of the electric field formulation (E-formulation) derived below that Mimp ∈ (L2 (Ω))3 . If the prescribed Mimp ∈ / (L2 (Ω))3 , we may still solve Maxwell’s equations with H(curl; Ω)-conforming finite elements for the magnetic field by using the H-formulation (2.3.2), or simply by prescribing an ˜ imp such that Mimp − M ˜ imp does not radiate outside the equivalent source M antenna [20]. Similarly, for the H-formulation, we shall assume that Jimp ∈ (L2 (Ω))3 . 2.3.1 E-Formulation First, we multiply the inverse of the magnetic permeability tensor by Fara¯ where F ∈ day’s law. Then, we multiply the resulting equation by ∇×F, HΓE (curl; Ω) = {F ∈ H(curl; Ω) : (n×F)|ΓE = 0 } is an arbitrary test function, and the symbol ’-’ denotes complex conjugate. Finally, integrating over the domain Ω, we arrive at the identity: Z

Ω

¯ ¯ −1 ∇×E) · (∇×F)dV (µ = −jω Z

Z

Ω

¯ dV H · (∇×F)

¯ dV . ¯ −1 Mimp ) · (∇×F) − (µ Ω

6

(2.9)

By integrating

Ω

obtain: Z

Ω

Z

¯ dV by parts, and by applying Ampere’s law, we H · (∇×F)

¯ dV = H · (∇×F) Z

Ω

Z

Ω

¯ dV − (∇×H) · F

¯ dV + ¯ + jω¯ (σ ǫ)E · F

Z

J

imp

Ω

Z

ΓN

¯ t dS = n ˆ×H · F

¯ dV − ·F

Z

ΓH

¯ t dS, n ˆ ×H · F

(2.10)

where Ft = F−(F· n ˆ )· n ˆ is the tangential component of vector F on ΓH , and n ˆ is the unit outward normal (with respect Ω if ΓH ⊂ ∂Ω) vector. Substitution of (2.10) into (2.9), together with equation (2.5) yields the following variational formulation:    Find E ∈ EΓE + HΓE (curl; Ω) such that:    Z Z    −1  ¯ ¯ dV = ¯  (µ ∇×E) · (∇×F) dV − (k¯2 E) · F    Ω Z Z Ω imp ¯ dV + jω Jimp ¯  ·F  −jω J ΓH · Ft dS     ΓH   Z Ω   −1  imp ¯ dV  ¯ M ) · (∇×F) − (µ ∀ F ∈ HΓE (curl; Ω)   

(2.11) ,

Ω

where k¯2 = ω 2¯ ǫ − jω σ ¯ , and EΓE is a lift (typically EΓE = 0) of the essential boundary condition data EΓE (denoted with the same symbol) related to Mimp ΓE by (2.7). Conversely, we can derive (2.1), (2.3), and (2.5) from the variational problem (2.11).

2.3.2 H-Formulation First, we multiply (σ ¯ + jω¯ ǫ)−1 by Ampere’s law. Then, we multiply the re¯ where F ∈ HΓ (curl; Ω) = {F ∈ H(curl; Ω) : sulting equation by ∇×F, H (n×F)|ΓH = 0 } is an arbitrary test function. Finally, by integrating over the domain Ω, we arrive at the identity: −jω

Z h

Ω

i

¯ (k¯2 )−1 (∇×H) · (∇×F)dV = −jω

Z h

Ω

Z

Ω

¯ dV E · (∇×F) i

¯ dV . (k¯2 )−1 Jimp · (∇×F)

7

(2.12)

By integrating

Ω

obtain: Z

Ω

Z

¯ dV by parts, and by applying Faraday’s law, we E · (∇×F)

¯ dV = E · (∇×F) −jω

Z

Ω

Z

Ω

¯ dV − (∇×E) · F

¯ dV − ¯ (µH) ·F

Z

imp

M

Ω

Z

ΓE

¯ t dS = n ˆ ×EΓE · F

¯ dV − ·F

Z

ΓE

¯ t dS. n ˆ×EΓE · F

(2.13)

Substitution of (2.13) into (2.12), together with equation (2.6) yields the following variational formulation:    Find H ∈ HΓH + HΓH (curl; Ω) such that:     Z h Z  i   −1  ¯ ¯ ¯ ¯ ¯  (σ + jω ǫ) (∇×H) · (∇×F)dV + jω (µH) · FdV    Ω Ω Z Z imp imp ¯ ¯ t dS+  − M · FdV + MΓE · F      ΓE   Z Ω i h    −1 imp ¯  ¯ + jω¯ · (∇×F)dV ∀ F ∈ HΓH (curl; Ω), (σ ǫ) J   

= (2.14)

Ω

where HΓH is a lift (typically HΓH = 0) of the essential boundary condition data HΓH (denoted with the same symbol) related to Jimp ΓH by (2.8).

2.4 Cylindrical Coordinates and Axisymmetric Problems

We consider a cylindrical coordinate system (ρ, φ, z). For the geophysical logging applications considered in this article, we assume that both the logging instrument and the rock formation properties are axisymmetric (invariant with respect to the azimuthal component φ) around the axis of the borehole. Under ˆ φ + zˆAz , this assumption, we obtain that for any vector field A = ρˆAρ + φA

∇×A = −ˆ ρ

1 ∂(ρAφ ) ∂Aφ ˆ ∂Aρ ∂Az + φ( − ) + zˆ . ∂z ∂z ∂ρ ρ ∂ρ 8

(2.15)

¯ and ¯ We also assume that tensors σ ¯ , µ, ǫ, describing the material properties of the medium, can be represented as:

¯= Ξ



 Ξρρ     0   



0 Ξρz 

Ξφφ 0

Ξzρ 0 Ξzz

      



=

 Ξρρ     0   



0 Ξρz  0 0

Ξzρ 0 Ξzz

|

{z

0       +0      



0 0

Ξφφ

   0   

¯ +Ξ ¯ , =Ξ ρ,z φ

(2.16)

0 0 0

}

¯ Ξ ρ,z



|

{z

¯ Ξ φ

}

¯ denotes either σ ¯ , µ, ¯ or ¯ where Ξ ǫ. According to this decomposition in terms ¯ , and Ξ ¯ ,respectively), we of the meridian and azimuthal “components” (Ξ ρ,z φ ¯ ρ,z , µ ¯ ρ,z , and ¯ define the meridian material properties σ ǫρ,z , and the azimuthal ¯ φ, µ ¯ φ , and ¯ ¯ ρ,z , material properties σ ǫφ . We also define k¯2 ρ,z = ω 2¯ ǫρ,z − jω σ ¯ φ. and k¯2 φ = ω 2¯ ǫφ − jω σ In order to simplify the notation, we introduce for every vector quantity F = (Fρ , Fφ , Fz ), its meridian part Fρ,z = (Fρ , 0, Fz ), and its azimuthal part Fφ = (0, Fφ , 0), so that F = Fρ,z ⊕ Fφ . 2.4.1 E-Formulation Next, we consider the space of all test functions F ∈ HΓE (curl; Ω) such that F = Fφ = (0, Fφ , 0). According to (2.15), ∇×F = −ˆ ρ

∂Fφ 1 ∂(ρFφ ) + zˆ . ∂z ρ ∂ρ

(2.17)

The variational formulation (2.11) reduces to a formulation in terms of the scalar field Eφ only, namely,   ˜ 1 (Ω) such that:  Find Eφ ∈ Eφ,ΓE + H  ΓE    Z Z    −1  ¯ ¯φ ¯  ( µ ∇×E ) · (∇× F ) dV − (k¯2 φ Eφ ) · F φ φ  ρ,z   Ω Ω Z Z imp imp ¯  −jω Jφ Fφ dV + jω Jφ,ΓH F¯φ dS      ΓH   Z Ω   −1  imp ¯ φ dV ˜ 1 (Ω) ,  ¯ ρ,z Mρ,z ) · F − (µ ∀ Fφ ∈ H  ΓE   Ω

9

dV = (2.18)

˜ 1 (Ω) = {Eφ : Eφ ∈ HΓ (curl; Ω)} = {Eφ ∈ L2 (Ω) : 1 Eφ + ∂Eφ ∈ where H ΓE E ρ ∂ρ ∂E φ ∈ L2 (Ω), Eφ |ΓE = 0}. Similarly, for a test function F = Fρ,z , the L2 (Ω) , ∂z variational problem (2.11) simplifies to:   ˜ Γ (curl; Ω) such that:   Find (Eρ , Ez ) ∈ EΓE + H E    Z Z    −1  ¯ ¯ ρ,z dV = ¯  (µφ ∇ × Eρ,z ) · (∇ × Fρ,z ) dV − (k¯2 ρ,z Eρ,z ) · F    Ω Ω Z Z imp ¯ imp ¯ imp ¯ imp ¯  −jω Jρ Fρ + Jz Fz dV + jω Jρ,Γ Fρ + Jz,Γ F dS   H H z    ΓH   Z Ω   −1 imp  ¯ ˜ Γ (curl; Ω) ,  ¯ − (µφ Mφ ) · Fρ,z dV ∀ (Fρ , Fz ) ∈ H  E  

(2.19)

Ω

˜ Γ (curl; Ω) = {(Eρ , Ez ) ∈ (L2 (Ω))2 : (∇×Eρ,z )|φ = ∂Eρ − ∂Ez ∈ where H E ∂z ∂ρ L2 (Ω) , (ˆ n×Eρ,z )|ΓE = 0}. In summary, problem (2.11) decouples into a system of two simpler problems described by (2.18) and (2.19).

¯ cannot be decomposed Remark 1. If |σφ,ρ |+|σρ,φ |+|σφ,z |+|σz,φ | = 6 0, then σ as the sum of its meridian and azimuthal “components” and, as a result, we cannot replace problem (2.11) with two simpler problems of lower dimension. Similarly, conditions |µφ,ρ| + |µρ,φ | + |µφ,z | + |µz,φ| = 0, and |ǫφ,ρ | + |ǫρ,φ | + |ǫφ,z | + |ǫz,φ | = 0 are essential in order to derive the formulations (2.18) and (2.19).

˜ Γ1 (Ω) It has been shown in [1] (Lemma 4.9) that the space H E 1 2 1 2 ˜ can also be expressed as HΓE (Ω) = {Eφ ∈ L (Ω) : Eφ ∈ L (Ω) , ∇(ρ,z) Eφ ∈ ρ 2 L (Ω), Eφ |ΓE = 0}.

Remark 2.

2.4.2 H-Formulation Using the same decomposition of test functions (i.e., F = Fφ , and F = Fρ,z ) for variational problem (2.14), we arrive at the following two decoupled 10

variational problems in terms of Hφ (2.20), and Hρ,z (2.21), respectively:   ˜ 1 (Ω) such that:  Find Hφ ∈ Hφ,ΓH + H  ΓH    Z h  i   −1  ¯ φ ) dV ¯ ρ,z + jω¯  ( σ ǫ ) ∇×H · (∇×F ρ,z φ    Ω Z Z Z imp ¯ imp ¯ ¯  ¯ Mφ,Γ Fφ dS  +jω (µφ Hφ ) · Fφ dV = − Mφ Fφ dV +  E    Ω Ω Γ  E  Z h i    −1 imp ¯ ˜ Γ1 (Ω)  ¯ ¯ ∀ Fφ ∈ H + (σ ρ,z + jω ǫρ,z ) Jρ,z · ∇×Fφ dV   H 

(2.20) ,

Ω

˜ 1 (Ω) = {Eφ : Eφ ∈ HΓ (curl; Ω)} = {Eφ ∈ L2 (Ω) : 1 Eφ + ∂Eφ ∈ where H ΓH H ρ ∂ρ ∂E φ L2 (Ω) , ∈ L2 (Ω), Eφ |ΓH = 0}. ∂z    ˜ Γ (curl; Ω) such that:  Find H = (Hρ , Hz ) ∈ HΓH + H  H   Z  h i    ¯ ρ,z ) dV ¯ φ + jω¯  (σ ǫφ )−1 ∇×Hρ,z · (∇ × F      Ω  Z Z    imp imp ¯  ¯  +jω (µρ,z Hρ,z ) · Fρ,z dV = − Mρ F¯ρ + Mz F¯z dV   Ω Z Ω  imp ¯ imp ¯   Mρ,ΓE Fρ + Mz,ΓE Fz dS +     ΓE  Z h  i   −1 imp  ¯ ρ,z ) dV ¯ φ + jω¯  + ( σ ǫ ) J · (∇ × F φ  φ     Ω     ˜ Γ (curl; Ω) ,  ∀ (Fρ , Fz ) ∈ H

(2.21)

H

˜ Γ (curl; Ω) = {(Eρ , Ez ) ∈ (L2 (Ω))2 , (∇×Eρ,z )|φ = ∂Eρ − ∂Ez ∈ where H H ∂z ∂ρ 2 L (Ω) , (ˆ n×E)|ΓH = 0}. From the formulation of problems (2.18) through (2.21), we remark the following: • Physically, the solution of problems (2.19), and (2.20) correspond to the T Eφ -mode (i.e. Eφ = 0), whereas the solution of problems (2.18), and (2.21) correspond to the T Mφ -mode (i.e. Hφ = 0). • The axis of symmetry is not a boundary of the original 3D problem, and therefore, a boundary condition (BC) is not needed to solve this problem. Nevertheless, formulations of problems (2.18) through (2.21) require the ˜ 1 (Ω), H ˜ 1 (Ω) and H ˜ Γ (curl; Ω), H ˜ Γ (curl; Ω) described use of spaces H ΓE ΓH E H above. The two former spaces involve the singular weight ρ1 , which implic11

itly requires a homogeneous Dirichlet boundary condition along the axis of symmetry. The two latter spaces are two-dimensional H(curl) spaces with no BC on the axis of symmetry. Thus, no BC is used to solve the problem. From the computational point of view, the effect of having no BC can be achieved by artificially adding a homogeneous natural (Neumann) BC.

3

Self-Adaptive Goal-Oriented hp-FEM

We are interested in solving the variational problems (2.11) and (2.14) (or alternatively, (2.18), (2.19), (2.20), and (2.21)), that we state here in terms of sesquilinear form b, and antilinear form f :     Find

E ∈ ED + V

   b(E, F)

(3.22)

= f (F) ∀F ∈ V ,

where, • • • •

ED is a lift of the essential (Dirichlet) BC. V is a Hilbert space. f ∈ V′ is an antilinear and continuous functional on V. b is a sesquilinear form. We have:

b(E, F) =

Z        Ω   |   Z         Ω   |

¯ dV − ¯ −1 ∇×E) · (∇×F) (µ {z

}

aE (E,F) −1

Z

Ω

|

¯ dV − (k¯2 ∇×E) · (∇×F) {z

}

aH (E,F)

¯ dV E-Form. k¯2 E · F {z

}

{z

}

cE (E,F)

Z

Ω

|

¯ dV H-Form. ¯ ·F µE cH (E,F)

, (3.23)

where sesquilinear forms aE , aH , cE and cH are Hermitian, continuous and semi-positive definite. We define an “energy” inner product on V as:

(E, F) :=

Z        Ω   |   Z         Ω   |

¯ dV + ¯ −1 ∇×E) · (∇×F) (µ {z

}

aE (E,F)

Z

Ω

|

¯ dV + (|k¯2 |−1 ∇×E) · (∇×F) {z

}

a ˜H (E,F)

¯ dV E-Form. |k¯2 |E · F Z

Ω

|

{z

c˜E (E,F)

}

¯ dV H-Form. ¯ ·F µE {z

cH (E,F)

}

, (3.24)

with the corresponding (energy) norm denoted by kEk. The symbol |k¯2 | 12

refers to a matrix of the same dimensions as k¯2 with its entries being the absolute value of the entries corresponding to k¯2 . At this point, we have assumed that all entries in the tensors included in eq. (3.24) are non-negative. 3.1 Representation of the Error in the Quantity of Interest Given an hp-FE subspace Vhp ⊂ V, we discretize (3.22) as follows:     Find

Ehp ∈ ED + Vhp

   b(Ehp , Fhp )

(3.25)

= f (Fhp ) ∀Fhp ∈ Vhp .

The objective of goal-oriented adaptivity 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 L, with L denoting a linear and continuous functional. By recalling the linearity of L, we have: Error of interest = L(E) − L(Ehp ) = L(E − Ehp ) = L(e) ,

(3.26)

where e = E−Ehp denotes the error function. By defining the residual rhp ∈ V′ as rhp (F) = f (F) − b(Ehp , F) = b(E − Ehp , F) = b(e, F), we seek the solution of the dual problem:     Find

¯ ∈V W

  ¯  b(F, W)

(3.27)

= L(F) ∀F ∈ V .

Problem (3.27) has a unique solution in V. The solution W, is usually referred to as the influence function. By discretizing (3.27) via, for example, Vhp ⊂ V, we obtain:     Find

¯ hp ∈ Vhp W

  ¯ hp )  b(Fhp , W

(3.28)

= L(Fhp ) ∀Fhp ∈ Vhp .

The definition of the dual problem plus the Galerkin orthogonality for the original problem imply the final representation formula for the error in the quantity of interest, namely, L(e) = b(e, W) = b(e, W − Fhp ) = ˜b(e, ǫ) . |

{z ǫ¯

}

13

At this point, Fhp ∈ Vhp is arbitrary, and ˜b(e, ǫ) = b(e, ǫ¯) denotes the bilinear form corresponding to the original sesquilinear form. Notice that, in practice, the dual problem is solved not for W but for its ¯ utilizing the bilinear form and not the sesquilinear form. complex conjugate W The linear system of equations is factorized only once, and the extra cost of solving (3.28) reduces to only one backward and one forward substitution (if a direct solver is used). Once the error in the quantity of interest has been determined in terms of the bilinear form ˜b, we wish to obtain a sharp upper bound for |L(e)| that depends upon the mesh parameters (element size h and order of approximation p) only locally. Then, a self-adaptive algorithm will be defined with the intent to minimize this bound. First, using a procedure similar to the one described in [7], we approximate E and W with fine grid functions E h , p+1 , W h , p+1 , which have been obtained 2 2 by solving the corresponding linear system of equations associated with the FE subspace V h , p+1 . In the remainder of this article, E and W will denote 2 the fine grid solutions of the direct and dual problems (E = E h , p+1 , and 2 W = W h , p+1, respectively), and we will restrict ourselves to discrete FE 2 spaces only. Next, we bound the error in the quantity of interest by a sum of element contributions. Let bK denote a contribution from element K to the sesquilinear form b. It then follows that |L(e)| = |b(e, ǫ)| ≤

X K

|bK (e, ǫ)| ,

(3.29)

where summation over K indicates summation over elements. Note again that e denotes now the difference between the fine and coarse grid solutions. Similarly ǫ is the difference between the fine grid approximation of the influence function and an arbitrary function Fhp ∈ Vhp . 3.2 Projection based interpolation operator Once we have a representation formula for the error in the quantity of interest in terms of the sum of element contributions given by (3.29), we wish to express this upper bound in terms of local quantities, i.e. in terms of quantities that do not vary globally when we modify the grid locally. For this purpose, we introduce the idea of projection-based interpolation operators. 14

First, in order to simplify the notation, we define the following three spaces of admissible solutions: • V = HD (curl; Ω), ˜ D (curl; Ω), and, • V2D = H ˜ 1 (Ω). • V 1D = H D 2D The corresponding hp-Finite Element spaces will be denoted by Vhp , Vhp , 1D and Vhp , respectively.

At this point, we introduce three projection-based interpolation operators that have been defined in [6,5], and used in [7,17] for the construction of the fully automatic energy-norm based hp-adaptive algorithm: • Πcurl,3D : V −→ Vhp , hp curl,2D 2D • Πhp : V2D −→ Vhp , and, 1D 1D 1D • Πhp : V −→ Vhp . We shall also consider three Galerkin projection operators: • Pcurl,3D : V −→ Vhp , hp curl,2D 2D , and, • Php : V2D −→ Vhp 1D 1D 1D • Php : V −→ Vhp . To further simplify the notation, we will utilize the unique symbol Πcurl hp to denote all projection based interpolation operators mentioned above. Depending upon the problem formulation (and corresponding space of admissible curl,3D for problems (2.11) and solutions), Πcurl hp should be understood as Πhp curl,2D 1D (2.14), Πhp for problems (2.19) and (2.21), or Πhp for problems (2.18) and curl,3D (2.20). Similarly, we will use the unique symbol Pcurl , hp to denote either Php curl,2D 1D Php , or Php . We denote Ehp = Pcurl hp E. Equation (3.29) then becomes |L(e)| ≤

X K

|bK (E, ǫ)| =

X K

curl curl |bK (E − Πcurl hp E, ǫ) + bK (Πhp E − Php E, ǫ)| .(3.30)

Given an element K, we conjecture that |bK (Πcurl hp E − Php E, ǫ)| will be negcurl ligible compared to |bK (E − Πhp E, ǫ)|. Under this assumption, we conclude that: |L(e)| 

X K

|bK (E − Πcurl hp E, ǫ)| . 15

(3.31)

In particular, for ǫ = W − Πcurl hp W, we have: |L(e)| 

X K

curl |bK (E − Πcurl hp E, W − Πhp W)| .

(3.32)

By applying Cauchy-Schwartz’s inequality, we obtain the next upper bound for |L(e)|: |L(e)| 

X K

k˜ ekK k˜ǫkK ,

(3.33)

˜ = E − Πcurl where e ǫ = W − Πcurl hp E, ˜ hp W, and k · kK denotes energy-norm k · k restricted to element K. 3.3 Fully Automatic Goal-Oriented hp-Refinement Algorithm We describe an hp self-adaptive algorithm that utilizes the main ideas of the fully automatic (energy-norm based) hp-adaptive algorithm described in [7,17]. We start by recalling the main objective of the self-adaptive (energy-norm based) hp-refinement strategy, which consists of solving the following maximization problem:    Find      ˜    hp

˜ an optimal hp-grid in the following sense:

= arg max

where

b hp

X K

2 curl 2 kE − Πcurl hp EkK − kE − Πhp b EkK

∆N

(3.34) ,

• E = E h , p+1 is the fine grid solution, and 2 • ∆N > 0 is the increment in the number of unknowns from grid hp to grid c hp.

Similarly, for goal-oriented hp-adaptivity, we propose the following algorithm based on estimate (3.33):     Find         ˜

˜ an optimal hp-grid in the following sense:

hp = arg max

          

b hp

X K

"

−

curl kE − Πcurl hp EkK · kW − Πhp WkK ∆N  curl curl kE − Πhp Ek · kW − Π Wk K K b b hp

∆N

16



(3.35) ,

where: • E = E h , p+1 and W = W h , p+1 are the fine grid solutions corresponding to 2 2 the direct and dual problems, and • ∆N > 0 is the increment in the number of unknowns from grid hp to grid c hp. The implementation of the goal-oriented hp-adaptive algorithm is based on the optimization procedure used for the energy-norm hp-adaptivity [7,17], which utilizes a multi-step approach (first optimization of edges, and then optimization of interior degrees of freedom). The subspace associated with an optimal FE grid is always contained in the subspace associated with the FE fine grid computed during the previous step.

3.4 Implementation details We outline shortly the main implementation details needed to extend the fully automatic (energy-norm based) hp-adaptive algorithm [7,17] to a fully automatic goal-oriented hp-adaptive algorithm. (1) First, the solution W of the dual problem on the fine grid is necessary. This goal can be attained either by using a direct (frontal) solver or an iterative (two-grid) solver (see [13]). (2) Subsequently, we treat both solutions as satisfying two different partial differential equations (PDE’s). We select the functions E and W as the solutions of the system of two PDE’s. (3) We proceed to redefine the evaluation of the error. The energy-norm error evaluation of a two-dimensional function is replaced with the product curl k E − Πcurl hp E k · k W − Πhp W k. (4) After these simple modifications, the energy-norm based self-adaptive algorithm may now be utilized as a self-adaptive goal-oriented hp algorithm.

4

Numerical Results

In this Section, we first describe a model problem corresponding to an alternating current (AC) through-casing resistivity tool (TCRT) in a borehole environment. Then, we utilize this problem to illustrate the main characteristics of the self-adaptive goal-oriented hp-FEM described in this article, namely: • The advantages of using goal-oriented hp-adaptivity as opposed to energy17

norm hp-adaptivity, • the advantages and disadvantages of using the azimuthal formulation (with continuous elements) as opposed to the meridian formulation (with edge elements), and • the final hp-grids, accuracy, and performance delivered by the self-adaptive goal-oriented hp-FEM.

4.1 Modeling of an AC Through-Casing Resistivity Tool (TCRT) in a Borehole Environment We consider a TCRT instrument composed of the following axisymmetric materials (all dimensions are given in cm): • one transmitter and two receiver coils defined on (1) ΩC1 = {(ρ, φ, z) : 0.08 m < ρ < 0.1 m , −0.01025 m < z < −0.00975 m}, (2) ΩC2 = {(ρ, φ, z) : 0.08 m < ρ < 0.1 m , 0.006375 m < z < 0.006625 m}, and, (3) ΩC3 = {(ρ, φ, z) : 0.08 m < ρ < 0.1 m , 0.008375 m < z < 0.008625 m}, respectively. The transmitter antenna is composed of a toroidal coil with an impressed volume magnetic current IM = 1/(πa2 ) Amperes, where a is the radius of the toroidal coil. The corresponding magnetic far-field solution in homogeneous and isotropic medium is given by (see [10]): −jkd j ρ ˆ + jωǫ)jk e H = φ(σ [1 − ] . 4πd kd d

(4.36)

This TCRT instrument moves along the vertical direction (z-axis) in a subsurface borehole environment composed of: • a borehole mud with resistivity 0.1 Ω·m defined on (1) ΩBH = {(ρ, φ, z) : ρ < 0.1 m} , • a casing, with resistivity 10−6 Ω·m defined on (1) ΩCS = {(ρ, φ, z) : 0.1 m ≤ ρ ≤ 0.1127 m} , and, • three formation materials of resistivities 100 Ω·m, 10000 Ω·m, and 1 Ω·m, defined on (1) ΩM1 = {(ρ, φ, z) : ρ ≥ 0.1127 m , (z < −0.5 m or z > 1 m)}, (2) ΩM2 = {(ρ, φ, z) : ρ ≥ 0.1127 m , −0.5 m ≤ z < 0 m}, and, (3) ΩM3 = {(ρ, φ, z) : ρ ≥ 0.1127 m , 0 m ≤ z ≤ 1 m}, respectively. Fig. 1 shows the geometry of the described logging instrument and borehole environment. 18

0.1m

0.1 Ohm−m 0.000001 Ohm−m

01 1010 1010 10 1010 1010 1010 1010 1010 1010 1010 1010 10 1010 1010 10 1010 10

0.15 m 0.20 m 0.65 m

0.5 m 0.5 m

z

100 Ohm−m

1 Ohm−m

10000 Ohm−m

100 Ohm−m

Fig. 1. 2D cross section of the geometry of a TCRT, composed of one transmitter and two receiver antennas (coils), a borehole, a 0.0127 m-thick steel casing, and four horizontal layers in the rock formation (with different resistivities).

From the engineering point of view, the main objective of these simulations is to compute the (normalized) difference of the vertical component of the electric field at the receiving coils. More precisely, we are interested in computing the quantity of interest Q given by  

Q(E) = 

1 |ΩC1 |

Z

ΩC1

Ez (l) dl −

1 |ΩC2 |

Z

ΩC2

 

Ez (l) dl /d(C1 , C2 ) ,

(4.37)

where d(C1 , C2 ) is the distance between the two receiver antennas. From the numerical point of view, solution of this TCRT problem is rather challenging for the following reasons: • It incorporates three singularities (of different strength) located at points where three different materials meet. • The resistivity between adjacent materials varies by as much as ten orders of magnitude. • The quantity of interest Q is typically several orders of magnitude smaller than the electric field measured at the receiver coils (see [14] for details). Thus, it is necessary to approximate the solution at the receiver coils very accurately, in order to obtain an acceptable error level in the quantity of interest. • Due to the presence of steel casing, we need to either consider a long computation domain in the vertical direction (3000 m or more) or to close the computational domain with an absorbing-type boundary condition (BC). Here, we consider a 4000 m domain in the vertical direction, and we impose 19

a homogeneous Dirichlet BC. 4.2 Numerical Comparison Between Energy-Norm and Goal-Oriented hpAdaptivity We consider the TCRT problem described in Section 4.1 operating at 10 Khz. In Fig. 2, we display the convergence history obtained using the self-adaptive energy-norm hp-algorithm [7]. The corresponding convergence history using the self-adaptive goal-oriented hp-algorithm described in this paper is shown in Fig. 3. From these figures, we make the following observations: • The approximate straight lines obtained in Figures 2 and 3 (in an algebraic vs. logarithmic scale) indicate the exponential convergence of both, the energy-norm hp-adaptivity (in terms of the relative energy-norm error vs. the number of unknowns) and the goal-oriented hp-adaptivity (in terms of upper bound (3.33) vs. the number of unknowns). The final slow-down in convergence associated with edge elements is explained in Section 4.3. • Exponential convergence in terms of the energy-norm does not imply an acceptable convergence in terms of the quantity of interest. Indeed, the final hp grid produced by the energy-norm adaptive algorithm exhibits a relative energy-norm error below 0.00001%, while the relative error in the quantity of interest remains large (over 200%). • The self-adaptive goal-oriented hp-FEM delivers exponential convergence rates in terms of estimate (3.33), which is an upper bound estimate for the error in the quantity of interest. Thus, the relative error in the quantity of interest rapidly approaches 0 (as the number of unknowns increases), even though the relative error in the quantity of interest may not decrease monotonically with the problem size. 4.3 Numerical Comparison Between Continuous and Edge Finite Elements At this point, we describe the main advantages and disadvantages of using edge elements with formulation (2.19) as opposed to continuous elements with formulation (2.20), for the case of simulating the TCRT problem considered in this paper: • The edge element formulation may entail a sudden slow-down in convergence due to round-off errors. As the frequency and/or finite element size decreases, the ratio between the L2 - and curl- contributions to the bilinear form b becomes smaller. In the limit, the L2 -contribution is not measurable with the computer (due to finite precision arithmetic), and the corresponding system of linear equations becomes singular (gradients become ill20

Continuous Elements

4

Edge Elements

4

10

10 Upper bound for |L(e)|/|L(u)| |L(e)|/|L(u)|

Energy−norm error |L(e)|/|L(u)|

2

2

10

Relative Error in %

Relative Error in %

10

0

10

−2

10

−4

10

−6

0

−2

10

−4

10

10

0

10

−6

1000 8000 27000 64000 1/3 Number of Unknowns N (scale N )

10

0

1000 8000 27000 64000 1/3 Number of Unknowns N (scale N )

Fig. 2. Convergence history for the TCRT problem operating at 10 Khz using energy-norm hp-adaptivity. The solid curves describe the relative energy-norm error, while dashed curves describe the relative error in the quantity of interest. Left panel: convergence history using formulation (2.20) with continuous (standard) elements. Right panel: convergence history using formulation (2.19) with edge elements. Frequency: 10 Khz

6

10

|L(e)|/|L(u)| (Cont. Elements) Upper bound for |L(e)|/|L(u)| (Cont. Elements) |L(e)|/|L(u)| (Edge Elements) Upper bound for |L(e)|/|L(u)| (Edge Elements)

4

Relative Error in %

10

2

10

0

10

−2

10

−4

10

1000

3375

8000 15625 27000 Number of Unknowns N (scale N1/3)

42875

Fig. 3. Convergence history for the TCRT problem operating at 10 Khz using goal-oriented hp-adaptivity. The solid curves describe the use of formulation (2.20) with continuous elements. The dashed curves describe the use of formulation (2.19) with edge elements. Light grey curves describe the relative error in the quantity of interest. Black curves describe the upper bound estimate (3.33).

21

defined). The slow-down in convergence observed in Figures 2 (right panel) and 3 correspond to the pre-limiting case, where the hp-grid contains elements for which the L2 -contribution is 14-15 orders of magnitude smaller than the corresponding curl-contribution. As we decrease the frequency, this slow-down in convergence occurs earlier, as shown in Fig. 4. The latter Frequency: 10 Hz

6

10

|L(e)|/|L(u)| (Cont. Elements) Upper bound for |L(e)|/|L(u)| (Cont. Elements) |L(e)|/|L(u)| (Edge Elements) Upper bound for |L(e)|/|L(u)| (Edge Elements)

4

Relative Error in %

10

2

10

0

10

−2

10

−4

10

1000

3375 8000 15625 27000 1/3 Number of Unknowns N (scale N )

42875

Fig. 4. Convergence history for the TCRT problem operating at 10 Hz using goal-oriented hp-adaptivity. Solid curves indicate the use of formulation (2.20) with continuous elements. Dashed curves indicate the use of formulation (2.19) with edge elements. Light grey curves display the relative error in the quantity of interest. Black curves display upper bound estimate (3.33).

stability problem can be overcome by explicitly re-enforcing the divergencefree condition (Gauss law) into the variational formulation at the expense of introducing a new variable (Lagrange multiplier). See [4] for details. • Continuous finite elements converge at a faster rate than edge elements (see Figures 3 and 4). This difference may be explained by the fact that an edge element of order p contains, for p > 1, more unknowns (2p2 + 2p) than a continuous element of order p (the number of unknowns is equal to p2 + 2p + 1). • For higher-order elements (p > 1), the number of unknowns associated with nodes located within the interior of an element is larger for edge elements than for continuous elements. Thus, for a fixed number of unknowns, a well-designed direct solver (performing static condensation on the interior of each element) should solve a problem using edge elements faster than using continuous elements 2 . 2

We do not display a convergence curve comparison against the CPU time due to

22

In summary, it is difficult to assess what type of elements (continuous or edge) are more suitable for simulating axisymmetric resistivity logging applications. In any case, by comparing results obtained from both the continuous and the edge element formulations we can verify the results, and perhaps utilize this comparison for constructing a-posteriori error estimates.

4.4 Final Grids, Accuracy, and Performance Delivered by the Self-Adaptive Goal-Oriented hp-FEM In the following, we display the final results obtained for the TCRT problem described in Section 4.1 operating at 10 Khz. Tables 1 and 2 display the final value of the quantity of interest Q. By using the fine grids as an error estimate for the coarse grids, we realize that the final coarse grids entail a relative error in the quantity of interest below 0.005%. By comparing fine grid results for continuous and edge elements, we estimate that both fine grids entail a relative error in the quantity of interest below 0.00002%. Table 1 Quantity of interest for continuous elements. Quantity of Interest

Real Part

Imag. Part

Final Coarse Grid

7.580492145E-09

-7.281872465e-09

Final Fine Grid

7.580470145E-09

-7.281954120e-09

Table 2 Quantity of interest for edge elements. Quantity of Interest

Real Part

Imag. Part

Final Coarse Grid

7.580304360e-09

-7.281686240e-09

Final Fine Grid

7.580469020e-09

-7.281954320e-09

The final hp-grids delivering a relative error in the quantity of interest below 0.1% for the TCRT problem operating at 10 Hz are displayed in Figures 5 (continuous elements) and 6 (edge elements). Table 3 displays the total CPU time associated with the self-adaptive goaloriented hp-FEM to solve the TCRT problem for 80 different vertical positions (spaced at 3 cm. intervals) of the logging instrument along the z-axis. For this performance analysis we use continuous elements and consider different antennas (a toroidal coil, a ring of vertical dipoles, a ring of horizontal dipoles, and a ring of electrodes) and different casing resistivities, ussing 10−6 Ω· m and 10−6 Ω· m, respectively. These results were obtained in an IBM Power 4 (1.3 Ghz) processor. The corresponding final logs are shown in Fig. 7. the lack of an efficient implementation for edge elements at this time.

23

f

f

Fig. 5. TCRT problem operating at 10Hz. Amplification (1.83m x 1.13m) of the final goal-oriented hp-grid in the vicinity of the two receiving antennas. Two circles indicate the location of the receiving antennas. This final goal-oriented hp-grid uses continuous elements and entails an error in the quantity of interest below 0.1%. Different colors indicate different polynomial orders of approximation, ranging from 1 to 8. Table 3 TCRT problem. Performance of the self-adaptive goal-oriented hp-FEM. 80 Vertical Positions

Electrical Resistivity of Casing 10−6 Ω · m

5

10−5 Ω · m

Toroid (10 Khz)

19’ 46”

16’ 28”

Ring of Vert. Dipoles (10 Khz)

22’ 47”

17’ 02”

Ring of Horiz. Dipoles (10 Khz)

19’ 25”

13’ 25”

Electrodes (0 Hz)

10’ 10”

08’ 35”

Summary and Conclusions

We have developed and succesfully tested a self-adaptive goal-oriented hp-FE algorithm that produces a sequence of hp-grids delivering exponential convergence rates in terms of a user-prescribed quantity of interest vs. the number 24

f

f

Fig. 6. TCRT problem operating at 10Hz. Amplification (1.83mx1.13m) of the final goal-oriented hp-grid in the vicinity of the two receiving antennas. Two circles indicate the location of the receiving antennas. This final goal-oriented hp-grid uses edge elements and entails an error in the quantity of interest below 0.1%. Different colors indicate different polynomial orders of approximation, ranging from 1 to 8.

of unknowns (as well as CPU time). The self-adaptive strategy combines the use of: • a projection based interpolation operator, which is independent of the material properties of the problem, and thus, it allows for a “simple” and problem-independent implementation, with • evaluation of element-based energy norms that incorporate material properties. This is essential for obtaining a sharp stable upper bound of the quantity of interest that we can use for minimization of the error. In the examples presented above, the ratio between the upper bound used for minimization and the quantity of interest typically ranges from 10 up to 100. Notice that if we had not included the material properties in the evaluation of the energy-norms, then this ratio would have become as large as 1011 , since the material properties vary by as much as ten orders of magnitude. Numerical results illustrate the performance and high-accuracy simulations obtained with the goal-oriented hp-FEM for problems that include several 25

Continuous Elements (10 Khz)

Position of the receiver antenna in the z−axis (in m)

2 Toroid Vert. dipoles (ring) Horiz. dipoles (ring) 1.5

1

0.5

0

−0.5

−1 −16 10

−14

−12

−10

−8

10 10 10 10 Amplitude of first difference of E (V/m2) z

Fig. 7. TCRT problem operating at 10KHz. Final log for different types of antennas. Results obtained using the self-adaptive goal-oriented hp-FEM.

singularities and large material contrasts. These results illustrate the necessity and importance of using goal-oriented algorithms for approximating an arbitrary quantity of interest of the solution. A comparison of continuous vs edge element discretizations for solving axisymmetric Maxwell’s equations reveals the advantages and disadvantages of each type of discretization. For the model problem considered in this paper, the use of edge elements provides a better initial approximation of the solution at low frequencies, but the rate of convergence is roughly half to that obtained with continuous elements. Since the influence function used by the self-adaptive goal-oriented hp-adaptive algorithm is approximated via Finite Elements, the numerical method presented in this article is problem independent, and therefore can be applied to FE discretizations of H 1 -, H(curl)-, and H(div)-spaces. 26

Acknowledgment This work was financially supported by Baker Atlas and The University of Texas at Austin’s Joint Industry Research Consortium on Formation Evaluation sponsored by Baker Atlas, BP, ConocoPhillips, ENI E&P, ExxonMobil, Halliburton Energy Services, Mexican Institute for Petroleum, Occidental Oil and Gas Corporation, Petrobras, Precision Energy Services, Schlumberger, Shell International E&P, Statoil, and Total. We would also like to acknowledge the expertise and technical advise received from L. Tabarovsky, A. Bespalov, T. Wang, and other members of the Science Department of Baker-Atlas.

References [1] F. Assous, C.(Jr.) Ciarlet, and S Labrunie, Theoretical tools to solve the axisymmetric Maxwell equations., Math. Meth. Appl. Sci., 25 (2002), pp. 49–78. [2] R. Becker and R. Rannacher, Weighted a posteriori error control in FE methods., in ENUMATH 97. Proceedings of the 2nd European conference on numerical mathematics and advanced applications held in Heidelberg, Germany, September 28-October 3, 1997. Including a selection of papers from the 1st conference (ENUMATH 95) held in Paris, France, September 1995. Singapore: World Scientific. 621-637 , Bock, Hans Georg (ed.) et al., 1998. [3] W. Cecot, W. Rachowicz, and L. Demkowicz, An hp-adaptive finite element method for electromagnetics. III: A three-dimensional infinite element for Maxwell’s equations., Int. J. Numer. Methods Eng., 57 (2003), pp. 899–921. [4] L. Demkowicz, hp-adaptive finite elements for time-harmonic Maxwell equations., in Topics in computational wave propagation. Direct and inverse problems. Berlin: Springer. Lect. Notes Comput. Sci. Eng. 31, 163-199 , Ainsworth, Mark (ed.) et al., 2003. [5]

, Finite element methods for Maxwell equations., Encyclopedia of Computational Mechanics, (eds. E. Stein, R. de Borst, T.J.R. Hughes), Wiley and Sons, 1 (2004).

[6] L. Demkowicz and A. Buffa, H 1 , H(curl), and H(div) conforming projection-based interpolation in three dimensions: quasi optimal p-interpolation estimates., Computer Methods in Applied Mechanics and Engineering, 194 (2005), pp. 267–296. [7] L. Demkowicz, W. Rachowicz, and Ph. Devloo, A fully automatic hpadaptivity., Journal of Scientific Computing, 17 (2002), pp. 117–142. [8] R. F. Harrington, Time-harmonic electromagnetic fields., McGraw-Hill, New York, 1961.

27

[9] V. Heuveline and R. Rannacher, Duality-based adaptivity in the hp-finite element method., Journal of Numerical Mathematics, 11 (2003), pp. 95–113. [10] J. R. Lovell, Finite element methods in resistivity logging, PhD thesis, Delft University of Technology, 1993. [11] J.T. Oden and S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method., Computers and Mathematics with Applications, 41 (2001), pp. 735–756. [12] M. Paraschivoiu and A. T. Patera, A hierarchical duality approach to bounds for the outputs of partial differential equations., Computer Methods in Applied Mechanics and Engineering, 158 (1998), pp. 389–407. [13] D. Pardo, Integration of hp-adaptivity with a two grid solver: applications to electromagnetics., PhD thesis, The University of Texas at Austin, April 2004. [14] D. Pardo, L. Demkowicz, and C. Torres-Verdin, A goal-oriented hp-adaptive finite element method with electromagnetic applications. Part I: electrostatics., ICES Report 04-57. To appear at Int. J. Numer. Methods Eng., 0 (2005), p. 0. [15] C. R. Paul and S. A. Nasar, Introduction to electromagnetic fields, McGrawHill, New York, 1982. [16] S. Prudhomme and J.T. Oden, On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors., Computer Methods in Applied Mechanics and Engineering, 176 (1999), pp. 313–331. [17] W. Rachowicz, D. Pardo, and L. Demkowicz, Fully automatic hpadaptivity in three dimensions., Tech. Report 04-22, ICES Report, 2004. [18] R. Rannacher and F.T. Suttmeier, A posteriori error control in finite element methods via duality techniques: application to perfect plasticity., Computational Mechanics, 21 (1998), pp. 123–133. [19] P. Solin and L. Demkowicz, Goal-oriented hp-adaptivity for elliptic problems., Computer Methods in Applied Mechanics and Engineering, 193 (2004), pp. 449–468. [20] J. Van Bladel, Singular electromagnetic fields and sources., Oxford University Press., New York, 1991.

28