A Appendix

# Robust A Posteriori Error Estimation for Finite Element Approximation to H(curl) Problem 1

## Abstract

In this paper, we introduce a novel a posteriori error estimator for the conforming finite element approximation to the problem with inhomogeneous media and with the right-hand side only in . The estimator is of the recovery type. Independent with the current approximation to the primary variable (the electric field), an auxiliary variable (the magnetizing field) is recovered in parallel by solving a similar problem. An alternate way of recovery is presented as well by localizing the error flux. The estimator is then defined as the sum of the modified element residual and the residual of the constitutive equation defining the auxiliary variable. It is proved that the estimator is approximately equal to the true error in the energy norm without the quasi-monotonicity assumption. Finally, we present numerical results for several interface problems.

## 1 Introduction

Let be a bounded and simply-connected polyhedral domain in with boundary and , and let be the outward unit vector normal to the boundary. Denote by the electric field, we consider the following model problem originated from a second order hyperbolic equation by eliminating the magnetic field in Maxwell’s equations:

 ⎧⎪ ⎪⎨⎪ ⎪⎩∇×(μ−1∇×u)+βu=f, in Ω,u×n=gD, on ΓD,(μ−1∇×u)×n=gN, % on ΓN, (1.1)

where is the curl operator; the , , and are given vector fields which are assumed to be well-defined on , , and , respectively; the is the magnetic permeability; and the depends on the electrical conductivity, the dielectric constant, and the time step size. Assume that the coefficients and are bounded below

 0<μ−10≤μ−1(x)and0<β0≤β(x)

for almost all .

The a posteriori error estimation for the conforming finite element approximation to the problem in (1.1) has been studied recently by several researchers. Several types of a posteriori error estimators have been introduced and analyzed. These include residual-based estimators and the corresponding convergence analysis (explicit [3, 10, 11, 12, 27, 29, 34], and implicit [18]), equilibrated estimators [4], and recovery-based estimators [6, 28]. There are four types of errors in the explicit residual-based estimator (see [3]). Two of them are standard, i.e., the element residual, and the interelement face jump induced by the discrepancy induced by integration by parts associated with the original equation in (1.1). The other two are also the element residual and the interelement face jump, but associated with the divergence of the original equation: , where is the divergence operator. These two quantities measure how good the approximation is in the kernel space of the curl operator.

Recently, the idea of the robust recovery estimator explored in [7, 8] for the diffusion interface problem has been extended to the interface problem in [6]. Instead of recovering two quantities in the continuous polynomial spaces like the extension of the popular Zienkiewicz-Zhu (ZZ) error estimator in [28], two quantities related to and are recovered in the respective - and -conforming finite element spaces. The resulting estimator consists of four terms similar to the residual estimator in the pioneering work [3] on this topic by Beck, Hiptmair, Hoppe, and Wohlmuth: two of them measure the face jumps of the tangential components and the normal component of the numerical approximations to and , respectively, and the other two are element residuals of the recovery type.

All existing a posteriori error estimators for the problem assume that the right-hand side is in or divergence free. This assumption does not hold in many applications (e.g. the implicit marching scheme mentioned in [19]). Moreover, two terms of the estimators are associated with the divergence of the original equation. In the proof, these two terms come to existence up after performing the integration by parts for the irrotational gradient part of the error, which lies in the kernel of the curl operator. One of the key technical tools, a Helmholtz decomposition, used in this proving mechanism, relies on being in , and fails if . In [12], the assumption that is weakened to being in the piecewise space with respect to the triangulation, at the same time, the divergence residual and norm jump are modified to incorporate this relaxation. Another drawback of using Helmholtz decomposition on the error is that it introduces the assumption of the coefficients’ quasi-monotonicity into the proof pipeline. An interpolant with a coefficient independent stability bound is impossible to construct in a “checkerboard” scenario (see [32] for diffusion case, and [6] for case). To gain certain robustness for the error estimator in the proof, one has to assume the coefficients distribution is quasi-monotone. However, in an earlier work of Chen, Xu, and Zou ([11]), it is shown that numerically this quasi-monotonicy assumption is more of an artifact introduced by the proof pipeline, at least for the irrotational vector fields. As a result, we conjecture that the divergence related terms should not be part of an estimator if it is appropriately constructed. In Section 5, some numerical justifications are presented to show the unnecessity of including the divergence related terms.

The pioneering work in using the dual problems for a posteriori error estimation dates back to [30]. In [30], Oden, Demkowicz, Rachowicz, and Westermann studied the a posteriori error estimation through duality for the diffusion-reaction problem. The finite element approximation to a dual problem is used to estimate the error for the original primal problem (diffusion-reaction). The result shares the same form to the Prague-Synge identity ([33]) for diffusion-reaction problem. The method presented in this paper may be viewed as an extension of the duality method in [30] to the interface problem. The auxiliary magnetizing field introduced in Section 3 is the dual variable resembling the flux variable in [30]. The connection is illustrated in details in Section 4.1.

Later, Repin ([31]) proposes a functional type a posteriori error estimator of problem, which can be viewed as an extension of the general approach in [30]. Repin et al ([26]) improve the estimate by assuming that the data is divergence free and the finite element approximation is in . In [31], the upper bound is established through integration by parts by introducing an auxiliary variable in an integral identity for . An auxiliary variable is recovered by globally solving an finite element approximation problem and is used in the error estimator. For the global lower bound, the error equation is solved globally in an conforming finite element space. Then the solution is inserted into the functional as the error estimator of which the maximizer corresponds to the solution to the error equation.

The purpose of this paper is to develop a novel a posteriori error estimator for the conforming finite element approximation to the problem in (2.1) that overcomes the above drawbacks of the existing estimators, e.g. the Helmholtz decomposition proof mechanism, restricted by the assumption that or divergence free, which brings in the divergence related terms. Specifically, the estimator studied in this paper is of the recovery type, requires the right-hand side merely having a regularity of , and has only two terms that measure the element residual and the tangential face jump of the original equation. Based on the current approximation to the primary variable (the electric field), an auxiliary variable (the magnetizing field) is recovered by approximating a similar auxiliary problem. To this end, a multigrid smoother is used to approximate this auxiliary problem, which is independent of the primary equation and is performed in parallel with the primary problem. The cost is the same order of complexity with computing the residual-based estimator, which is much less than solving the original problem.

An alternate route is illustrated as well in Section 3.2 by approximating a localized auxiliary problem. While embracing the locality, the parallel nature using the multigrid smoother is gone. The recovery through approximating localized problem requires the user to provide element residual and tangential face jump of the numerical magnetizing field based on the finite element solution of the primary equation. The estimator is then defined as the sum of the modified element residual and the residual of the auxiliary constitutive equation. It is proved that the estimator is equal to the true error in the energy norm globally. Moreover, in contrast to the mechanism of the proof using Helmholtz decomposition mentioned previously, the decomposition is avoided by using the joint energy norm. As a result, the new estimator’s reliability does not rely on the coefficients distribution (Theorem 4.2).

Meanwhile, in this paper, the method and analysis extend the functional-type error estimator in [31] to a more pragmatic context by including the mixed boundary conditions, and furthermore, the auxiliary variable is approximated by a fast multigrid smoother, or by solving a localized problem on vertex patches, to avoid solving a global finite element approximation problem.

Lastly, in order to compare the new estimator introduced in this paper with existing estimators, we present numerical results for intersecting interface problems. When , the mesh generated by our indicator is much more efficient than those by existing indicators (Section 5).

## 2 Primal Problem and The Finite Element Approximation

Denote by the space of the square integrable vector fields in equipped with the standard norm: , where denotes the standard inner product over an open subset , when , the subscript is dropped for and . Let

 H(curl;Ω):={v∈L2(Ω):∇×v∈L2(Ω)},

which is a Hilbert space equipped with the norm

 ∥v∥H(curl)=(∥v∥2+∥∇×v∥2)1/2.

Denote its subspaces by

 HB(curl;Ω):={v∈H(curl;Ω):v×n=gB on ΓB}and\vboxto−0.43pt∘\vssHB(curl;Ω):={v∈HB(curl;Ω):gB=0}

for or .

For any , multiplying the first equation in (1.1) by a suitable test function with vanishing tangential part on , integrating over the domain , and using integration by parts formula for -regular vector fields (e.g. see [2]), we have

 (f,v) = (∇×(μ−1∇×u),v)+(βu,v) = (μ−1∇×u,∇×v)+(βu,v)−∫ΓNgN⋅vdS.

Then the weak form associated to problem (1.1) is to find such that

 Aμ,β(u,v)=fN(v),∀v∈\vboxto−0.43pt∘\vssHD(curl;Ω), (2.1)

where the bilinear and linear forms are given by

 Aμ,β(u,v)=(μ−1∇×u,∇×v)+(βu,v)andfN(v)=(f,v)+⟨gN,v⟩ΓN,

respectively. Here, denotes the duality pair over . Denote by

 Missing dimension or its units for \mskip

the “energy” norm induced by the bilinear form .

###### Theorem 2.1.

Assume that , , and . Then the weak formulation of (1.1) has a unique solution satisfying the following a priori estimate

 Missing dimension or its units for \mskip (2.2)
###### Proof.

For the notations and proof, see the Appendix A. ∎

### 2.1 Finite Element Approximation

For simplicity of the presentation, only the tetrahedral elements are considered. Let be a finite element partition of the domain . Denote by the diameter of the element . Assume that the triangulation is regular and quasi-uniform.

Let where is the space of polynomials of degree less than or equal to . Let and be the spaces of homogeneous polynomials of scalar functions and vector fields. Denote by the first or second kind Nédélec elements (e.g. see [24, 25])

 Missing dimension or its units for \hskip

for , respectively, where the local Nédélec elements are given by

 NDk,1(K)={p+s:p∈Pk(K),s∈˜Pk+1(K) such that s⋅x=0} and NDk,2(K)={p+∇s:p∈NDk,1(K),s∈˜Pk+2(K)}.

For simplicity of the presentation, we assume that both boundary data and are piecewise polynomials, and the polynomial extension (see [14]) of the Dirichlet boundary data as the tangential trace is in . Now, the conforming finite element approximation to (1.1) is to find such that

 Aμ,β(uT,v)=fN(v),∀v∈NDk∩\vboxto−0.43pt∘\vssHD(curl;Ω). (2.3)

Assume that and are the solutions of the problems in (1.1) and (2.3), respectively, and that , (When the regularity assumption is not met, one can construct a curl-preserving mollification, see [16]), by the interpolation result from [24] Chapter 5 and Céa’s lemma, one has the following a priori error estimation:

 |||u−uT|||μ,β≤Chk+1(∥u∥Hk+1(Ω)+∥∇×u∥Hk+1(Ω)), (2.4)

where is a positive constant independent of the mesh size .

## 3 Auxiliary Problem of Magnetizing Field

### 3.1 Recovery of the magnetizing field

Introducing the magnetizing field

 σ=μ−1∇×u, (3.1)

then the first equation in (1.1) becomes

 ∇×σ+βu=f, in Ω. (3.2)

The boundary condition on may be rewritten as follows

 σ×n=gN, on ΓN.

For any , multiplying equation (3.2) by , integrating over the domain , and using integration by parts and (3.1), we have

 (β−1f,∇×τ) = (β−1∇×σ,∇×τ)+(u,∇×τ) = (β−1∇×σ,∇×τ)+(∇×u,τ) +∫ΓD(u×n)⋅τds−∫ΓNu⋅(τ×n)ds = (β−1∇×σ,∇×τ)+(μσ,τ)+∫ΓDgD⋅τds.

Hence, the variational formulation for the magnetizing field is to find such that

 Aβ,μ(σ,τ)=fD(τ),∀τ∈\vboxto−0.43pt∘\vssHN(curl;Ω), (3.3)

where the bilinear and linear forms are given by

 Aβ,μ(σ,τ)=(β−1∇×σ,∇×τ)+(μσ,τ)andfD(τ)=(β−1f,∇×τ)−⟨gD,τ⟩ΓD,

respectively. The natural boundary condition for the primary problem becomes the essential boundary condition for the auxiliary problem, while the essential boundary condition for the primary problem is now incorporated into the right-hand side and becomes the natural boundary condition. Denote the “energy” norm induced by by

 |||τ|||β,μ=√Aβ,μ(τ,τ).
###### Theorem 3.1.

Assume that , , and . Then problem (3.3) has a unique solution satisfying the following a priori estimate

 |||σ|||β,μ≤∥β−1/2f∥+∥∥gD∥∥1/2,μ,β,ΓD+∥∥gN∥∥−1/2,β,μ,ΓN. (3.4)
###### Proof.

The theorem may be proved in a similar fashion as Theorem 2.1. ∎

Similarly to that for the essential boundary condition, it is assumed that the polynomial extension of the Neumman boundary data as the tangential trace is in as well. Now, the conforming finite element approximation to (3.3) is to find such that

 Aβ,μ(σT,τ)=fD(τ),∀τ∈NDk∩\vboxto−0.43pt∘\vssHN(curl;Ω). (3.5)

Assume that and are the solutions of the problems in (3.1) and (3.5), respectively, and that , , one has the following a priori error estimation similar to (2.4)

 Missing or unrecognized delimiter for \Big (3.6)

The a priori estimate shows that heuristically, for the auxiliary magnetizing field , using the same order -conforming finite element approximation spaces with the primary variable may be served as the building blocks for the a posteriori error estimation.

### 3.2 Localization of the recovering procedure

The localization of the recovery of for this new recovery shares similar methodology with the one used in the equilibrated flux recovery (see [4, 5]). However, due to the presence of the -term, exact equilibration is impossible due to several discrepancies: if and are in Nédélec spaces of the same order; If is used for and for , the inter-element continuity conditions come into the context in that , which has different inter-element continuity requirement than . Due to these two concerns, the local problem is approximated using a constraint -minimization.

Let be the correction from to the true magnetizing field: . Now can be decomposed using a partition of unity: let be the linear Lagrange nodal basis function associated with a vertex , which is the collection of all the vertices,

 σΔ=∑z∈NhσΔz, with σΔz:=λzσΔ. (3.7)

Denote . Let the vertex patch , where is the collection of vertices of element . Then the following local problem is what the localized magnetizing field correction satisfies:

 {μσΔz−∇×ez=−∇λz×e, in K⊂ωz,∇×σΔz+βez=λzrK+∇λz×(μ−1∇×e), in K⊂ωz, (3.8)

with the following jump condition on each interior face , and boundary face :

 {[[σΔz×nF]]\raisebox−2.0pt$F$=−λzjt,F, on F∈Fz,σΔz×nF=0, on F⊂∂ωz. (3.9)

The element residual is , and the tangential jump is .

To find the correction, following piecewise polynomial spaces are defined:

 NDk−1(ωz)={τ∈L2(ωz):τ∣∣\raisebox−0.5pt$K$∈NDk(K),∀K⊂ωz}, (3.10) Wk(Fz)={τ∈L2(Fz):τ∣∣\raisebox−0.5pt$F$∈RTk(F),∀F∈Fz; τ∣∣\raisebox−0.5pt$Fi$⋅(tij×ni)=τ∣∣\raisebox−0.5pt$Fj$⋅(tij×nj),∀Fi,Fj∈Fz,∂Fi∩∂Fj=eij}, Hz={τ∈NDk−1(ωz):[[τ×nF]]\raisebox−2.0pt$F$=−¯jF,z∀F∈Fz}, and H0,z={τ∈Hz:τ×nF∣∣\raisebox−0.5pt$F$=0,∀F⊂ωz}.

Here is the planar Raviart-Thomas space on a given face , of which the degrees of freedom can be defined using conormal of an edge with respect to the face normal . For example, is the unit tangential vector of edge joining face and , then the conormal vector of with respect to face is . can be viewed as the trace space of the broken Nédélec space . For detail please refer to Section 4 and 5 in [13].

To approximate the local correction for magnetizing field, and are projected onto proper piecewise polynomial spaces. To this end, let

 ¯¯¯rK,z:=∏K(λzrK),%and¯jF,z:=∏F(λzjt,F), (3.11)

where is the projection onto the space , and is the projection onto the space . Dropping the uncomputable terms in (3.8), and using (3.9) as a constraint, the following local -minimization problem is to be approximated:

 Missing \left or extra \right (3.12)

The hybridized problem associated with above minimization is obtained by taking variation with respect to of the functional by the tangential face jump as a Lagrange multiplier:

 J∗z(σΔz,\raisebox−0.5pt$T$,ξ) :=12∥∥μ1/2σΔz,\raisebox−0.5pt$T$−μ−1/2∇×ez∥∥2ωz (3.13) +12∥∥β−1/2(∇×σΔz,\raisebox−0.5pt$T$+βez−¯¯¯rK,z)∥∥2ωz +∑F∈Fz([[σΔz,\raisebox−0.5pt$T$×nF]]\raisebox−2.0pt$F$+¯jF,z,ξ)F.

For any , using the fact that , and on

 0= (μ1/2σΔz,\raisebox−0.5pt$T$−μ−1/2∇×ez,μ1/2τ)ωz (3.14) +(β−1/2∇×σΔz,\raisebox−0.5pt$T$+β1/2ez−β−1/2¯¯¯rK,z,β−1/2∇×τ)ωz +∑F∈Fz([[τ×nF]]\raisebox−2.0pt$F$,ξ)F = Missing or unrecognized delimiter for \bigr +∑F∈Fz([[τ×nF]]\raisebox−2.0pt$F$,ξ−ez)F−(β−1¯¯¯rK,z,∇×τ)ωz.

As a result, the local approximation problem is:

 Missing or unrecognized delimiter for \big (3.15)

wherein the local bilinear forms are defined as follows:

 Aβ,μ;z(σ,τ) :=(β−1∇×σ,∇×τ)ωz+(μσ,τ)ωz, (3.16) and Bz(τ,γ) :=∑F∈Fz([[τ×nF]]\raisebox−2.0pt$F$,γ)F.
###### Proposition 3.2.

Problem (3.15) has a unique solution.

###### Proof.

For a finite dimensional problem, uniqueness implies existence. It suffices to show that letting both the right hand sides be zeros results trivial solution. First by for any (direct implication of Proposition 4.3 and Theorem 4.4 in [13]), setting in the second equation of (3.15) immediately implies that . As a result, . Now let in the first equation of (3.15), since induces a norm in , . For , it suffices to show that on each if

 Extra open brace or missing close brace

Using Theorem 4.4 in [13], if is non-trivial and satisfies above equation, there always exists a such that . As a result, , which is a contradiction. Thus, the local problem (3.15) is uniquely solvable. ∎

With the local correction to the magnetizing field, for all , computed above, let

 ˜σΔK,T=∑z∈N(K)σΔz,\raisebox−0.5pt$T$,and˜σΔT=∑z∈NσΔz,\raisebox−0.5pt$T$, (3.17)

then the recovered magnetizing field is

 ˜σT=˜σΔT+μ−1∇×uT. (3.18)

## 4 A Posteriori Error Estimator

In this section, we study the following a posteriori error estimator:

 η=(∑K∈Tη2K)1/2

where the local indicator is defined by

 ηK=(∥∥μ−1/2(μσT−∇×uT)∥∥2K+∥∥β−1/2(∇×σT+βuT−f)∥∥2K)1/2. (4.1)

It is easy to see that

 η=(∥∥μ−1/2(μσT−∇×uT)∥∥2+∥∥β−1/2(∇×σT+βuT−f)∥∥2)1/2. (4.2)

The and are the finite element approximations in problems (2.3) and (3.5) respectively.

With the locally recovered , the local error indicator and the global error estimator are defined in the same way as (4.1) and (4.2):

 ˜ηK=(∥∥μ−1/2(μ˜σT−∇×uT)∥∥2K+∥∥β−1/2(∇×˜σT+βuT−f)∥∥2K)1/2, (4.3)

and

 ˜η=(∥∥μ−1/2(μ˜σT−∇×uT)∥∥2+∥∥β−1/2(∇×˜σT+βuT−f)∥∥2)1/2. (4.4)
###### Remark 4.1.

In practice, does not have to be the finite element solution of a global problem. In the numerical computation, the Hiptmair-Xu multigrid preconditioner in [20] is used for discrete problem (3.5) with two multigrid V-cycles for each component of the vector Laplacian, and two multigrid V-cycles for the kernel part of the curl operator. The used to evaluate the estimator is the PCG iterate. The computational cost is the same order with computing the explicit residual based estimator in [3].

Generally speaking, to approximate the auxiliary problem, the same black-box solver for the original problem can be applied requiring minimum modifications. For example, if the BoomerAMG in hypre ([17, 22]) is used for the discretizations of the primary problem, then the user has to provide exactly the same discrete gradient matrix and vertex coordinates of the mesh, and in constructing the the HX preconditioner, the assembling routines for the vector Laplacian and scalar Laplacian matrices can be called twice with only the coefficients input switched.

###### Theorem 4.2.

Locally, the indicator and both have the following efficiency bound

 Missing dimension or its units for \mskip (4.5)

for all . The estimator and satisfy the following global upper bound

 Missing dimension or its units for \mskip (4.6)
###### Proof.

Denote the true errors in the electric and magnetizing fields by

 e=u−uT,andE=σ−σT,

respectively. It follows from (3.1), (3.2), and the triangle inequality that

 η2K = ∥∥μ1/2E−μ−1/2∇×e∥∥2K+∥∥β−1/2∇×E+β1/2e∥∥2K (4.7) ≤ (∥∥μ1/2E∥∥2K+∥∥μ−1/2∇×e∥∥2K+∥∥β−1/2∇×E∥∥2K+∥∥β1/2e∥∥2K) = (|||e|||2μ,β,K+|||E|||2β,μ,K),

which implies the validity of (4.5) for . For , the exact same argument follows except by switching by locally recovered . To prove the global identity in (4.6), summing (4.7) over all gives

 η2 = ∥∥μ1/2E−μ−1/2∇×e∥∥2+∥∥β−1/2∇×E+β1/2e∥∥2 = Missing dimension or its units for \mskip

Now, (4.6) follows from the fact that

 −(E,∇×e)+(∇×E,e)=0.

Lastly, the global upper bound for the locally recovered follows from the fact that