# An equilibrated a posteriori error estimator for arbitrary-order Nédélec elements for magnetostatic problems

We present a novel a posteriori error estimator for Nédélec elements for magnetostatic problems that is constant-free, i.e. it provides an upper bound on the error that does not involve a generic constant. The estimator is based on equilibration of the magnetic field and only involves small local problems that can be solved in parallel. Such an error estimator is already available for the lowest-degree Nédélec element [D. Braess, J. Schöberl, Equilibrated residual error estimator for edge elements, Math. Comp. 77 (2008)] and requires solving local problems on vertex patches. The novelty of our estimator is that it can be applied to Nédélec elements of arbitrary degree. Furthermore, our estimator does not require solving problems on vertex patches, but instead requires solving problems on only single elements, single faces, and very small sets of nodes. We prove reliability and efficiency of the estimator and present several numerical examples that confirm this.

## Authors

• 3 publications
• 5 publications
• 11 publications
04/17/2020

### A polynomial-degree-robust a posteriori error estimator for Nédélec discretizations of magnetostatic problems

We present an equilibration-based a posteriori error estimator for Nédél...
09/29/2020

### Residual-based a posteriori error estimates for 𝐡𝐩-discontinuous Galerkin discretisations of the biharmonic problem

We introduce a residual-based a posteriori error estimator for a novel h...
04/23/2020

### Adaptive virtual element methods with equilibrated flux

We present an hp-adaptive virtual element method (VEM) based on the hype...
05/06/2020

### A Residual Based A Posteriori Error Estimators for AFC Schemes for Convection-Diffusion Equations

In this work, we propose a residual-based a posteriori error estimator f...
04/20/2020

### Residual-type a posteriori error analysis of HDG methods for Neumann boundary control problems

We study a posteriori error analysis of linear-quadratic boundary contro...
06/17/2021

### A posteriori estimator for the adaptive solution of a quasi-static fracture phase-field model with irreversibility constraints

Within this article, we develop a residual type a posteriori error estim...
10/19/2020

### A posteriori error estimates by weakly symmetric stress reconstruction for the Biot problem

A posteriori error estimates are constructed for the three-field variati...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

We consider an a posteriori error estimator for finite element methods for solving equations of the form . These equations are related to magnetostatics but also appear in eddy current models for non-conductive media.

The first a posteriori error estimator in this context was introduced and analysed in [4]. It is a residual-type estimator and provides bounds of the form

 c0estimator≤error ≤c1estimator,

up to some higher-order data oscillation terms, where are positive constants that do not depend on the mesh resolution. Similar bounds can be obtained by hierarchical error estimators; see, e.g., [5], under the assumption of a saturation condition, and by Zienkiewicz–Zhu-type error estimators; see, e.g., [17]. A drawback of these estimators is that the constants and are usually unknown, resulting in significant overestimation or underestimation of the real error.

Equilibration-based error estimators can circumvent this problem. Often attributed to Prager and Synge [18], these estimators have become a major research topic; for a recent overview, see, for example, [13] and the references therein. An equilibration-based error estimator was introduced for magnetostatics in [7] and provides bounds of the form

 c0estimator≤error ≤estimator

up to some higher-order data oscillation terms. In other words, it provides a constant-free upper bound on the error. A different equilibration-based error estimator for magnetostatics was introduced in [19] and, for an eddy current problem, in [11, 10]. Constant-free upper bounds are also obtained by the functional estimate in [16], when selecting a proper function in their estimator, and by the recovery-type error estimator in [8], in case the equations contain an additional term , with .

A drawback of the estimators in [19, 11, 10, 16] is that they require solving a global problem. The estimator in [7], on the other hand, only involves solving local problems related to vertex patches. However, the latter estimator is defined for Nédélec elements of the lowest degree only. In this paper, we present a new equilibration-based constant-free error estimator that can be applied to Nédélec elements of arbitrary degree. Furthermore, our estimator involves solving problems on only single elements, single faces, and very small sets of nodes.

The paper is constructed as follows: We firstly introduce a finite element method for solving magnetostatic problems in Section 2. We then derive our error estimator step by step in Section 3, with a summary given in Section 3.3, and prove its reliability and efficiency in Section 3.4. Numerical examples confirming the reliability and efficiency of our estimator are presented in Section 4, and an overall summary is given in Section 5.

## 2. A finite element method for magnetostatic problems

Let be an open, bounded, simply connected, polyhedral domain with a connected Lipschitz boundary . In case of a linear, isotropic medium and a perfectly conducting boundary, the static magnetic field satisfies the equations

 ∇×H =j in Ω, ∇⋅μH =0 in Ω, ^n⋅μH =0 on ∂Ω,

where

is the vector of differential operators

, and denote the outer- and inner product, respectively, (therefore, and are the curl- and divergence operator, respectively), denotes the outward pointing unit normal vector, , with for some positive constants and , is a scalar magnetic permeability, and is a given divergence-free current density. The first equality is known as Ampères law and the second as Gauss’s law for magnetism.

These equations can be solved by writing , where is a vector potential, and by solving the following problem for :

 (1a) ∇×(μ−1∇×u) =j in Ω, (1b) ∇⋅u =0 in Ω, (1c) ^n×u =0 on ∂Ω.

The second condition is only added to ensure uniqueness of and is known as Coulomb’s gauge.

Now, for any domain , let denote the standard Lebesque space of square-integrable vector-valued functions equipped with norm and inner product , and define the following Sobolev spaces:

 H1(Ω) :={ϕ∈L2(Ω)|∇ϕ∈L2(Ω)3}, H10(Ω) :={ϕ∈H1(Ω)|ϕ=0 on ∂Ω}, H(curl;Ω) :={u∈L2(Ω)3|∇×u∈L2(Ω)3}, H0(curl;Ω) :={u∈H(curl;Ω)|^n×u=0 on ∂Ω}, H(div;Ω) :={u∈L2(Ω)3|∇⋅u∈L2(Ω)}, H(div0;Ω) :={u∈H(div;Ω)|∇⋅u=0}.

The weak formulation of problem (1) is finding such that

 (2) (μ−1∇×u,∇×w)Ω =(j,w)Ω ∀w∈H0(curl;Ω),

which is a well-posed problem [15, Theorem 5.9].

The solution of the weak formulation can be approximated using a finite element method. Let be a tetrahedron and define to be the space of polynomials on of degree or less. Also, define the Nédélec space of the first kind and the Raviart-Thomas space by

 Rk(T) Dk(T)

Finally, let denote a tessellation of into tetrahedra with a diameter smaller than or equal to , let , , and denote the discontinuous spaces given by

 P−1k(Th) :={ϕ∈L2(Ω)|ϕ|T∈Pk(T) for all% T∈Th}, R−1k(Th) :={u∈L2(Ω)3|u|T∈Rk(T) for all T∈Th}, D−1k(Th) :={u∈L2(Ω)3|u|T∈Dk(T) for all T∈Th},

and define

 Pk(Th) :=P−1k(Th)∩H1(Ω), Pk,0(Th) :=P−1k(Th)∩H10(Ω), Rk(Th) :=R−1k(Th)∩H(curl;Ω), Rk,0(Th) :=R−1k(Th)∩H0(curl;Ω), Dk(Th) :=D−1k(Th)∩H(div;Ω).

We define the finite element approximation for the magnetic vector potential as the vector field that solves

 (3a) (μ−1∇×uh,∇×w)Ω =(j,w)Ω ∀w∈Rk,0(Th), (3b) (uh,∇ψ)Ω =0, ∀ψ∈Pk,0(Th).

The approximation of the magnetic field is then given by

 Hh:=μ−1∇×uh,

which converges quasi-optimally as the mesh width tends to zero [15, Theorem 5.10].

In the next section, we show how we can obtain a reliable and efficient estimator for .

## 3. An equilibration-based a posteriori error estimator

We follow [7] and present an a posteriori error estimator that is based on the following result.

###### Theorem 3.1 ([7, Thm. 10]).

Let be the solution to (2), let be the solution of (3), and set and . If satisfies the equilibrium condition

 (4) ∇×~H =j,

then

 (5) ∥μ1/2(H−Hh)∥Ω ≤∥μ1/2(~H−Hh)∥Ω.
###### Proof.

The result follows from the orthogonality of and :

 =(μ1/2(~H−H),μ−1/2∇×(u−uh))Ω =(~H−H,∇×(u−uh))Ω =(∇×(~H−H),u−uh)Ω =(j−j,u−uh)Ω =0

and Pythagoras’s theorem

 (6) ∥μ1/2(~H−Hh)∥2Ω =∥μ1/2(~H−H)∥2Ω+∥μ1/2(H−Hh)∥2Ω.

###### Remark 3.2.

Equation (6) is also known as a Prager–Synge type equation and obtaining an error estimator from such an equation is also known as the hypercircle method. Furthermore, equation (4) is known as the equilibrium condition and using the numerical approximation to obtain a solution to this equation is called equilibration of .

###### Corollary 3.3.

Let be the solution to (2), let be the solution of (3), set and , and let be the discrete current distribution. If satisfies the (residual) equilibrium condition

 (7) ∇×~HΔ =j−jh,

which is an identity of distributions, then

 (8) ∥μ1/2(H−Hh)∥Ω ≤∥μ1/2~HΔ∥Ω.
###### Proof.

Since (7) is an identity of distributions, we can equivalently write

 ⟨∇×~HΔ,w⟩ =⟨j−jh,w⟩ ∀w∈C∞0(Ω)3,

where denotes the application of a distribution to a function in . Now, set . Using the definition , we obtain

 ⟨∇×~H,w⟩ =⟨∇×~HΔ+∇×Hh,w⟩ =⟨j−jh+∇×Hh,w⟩ =⟨j,w⟩ ∀w∈C∞0(Ω)3.

From this, it follows that , so is in and satisfies equilibrium condition (4). Inequality (8) then follows from Theorem 3.1. ∎

From Corollary 3.3, it follows that a constant-free upper bound on the error can be obtained from any field that satisfies (7).

An error estimator of this type was first introduced in [7], where it is referred to as an equilibrated residual error estimator. There, is decomposed into a sum of local divergence-free current distributions that have support on only a single vertex patch. The error estimator is then obtained by solving local problems of the form for each vertex patch and by then taking the sum of all local fields . It is, however, not straightforward to decompose into local divergence-free current distributions. An explicit expression for is given in [7] for the lowest-degree Nédélec element, but this expression cannot be readily extended to basis functions of arbitrary degree.

Here, we instead present an error estimator based on equilibration condition (7) that can be applied to elements of arbitrary degree. Furthermore, instead of solving local problems on vertex patches, our estimator requires solving problems on only single elements, single faces, and small sets of nodes. The assumptions and a step-by-step derivation of the estimator are given in Sections 3.1 and 3.2 below, a brief summary is given in Section 3.3, and reliability and efficiency are proven in Section 3.4.

### 3.1. Assumptions

In order to compute the error estimator, we use polynomial function spaces of degree , where denotes the degree of the finite element approximation , and assume that:

• The magnetic permeability is piecewise constant. In particular, the domain can be partitioned into a finite set of polyhedral subdomains such that is constant on each subdomain . Furthermore, the mesh is assumed to be aligned with this partition so that is constant within each element.

• The current density is in .

Although assumption A2 does not hold in general, we can always replace by a suitable projection by taking, for example,

as the standard quasi-interpolation operator of the

space. The error is in that case bounded by

 ∥μ1/2(H−Hh)∥2Ω ≤∥μ1/2~HΔ∥2Ω+C∥j−πhj∥Ω,

where is some positive constant that does not depend on the mesh width . If is sufficiently smooth, then the data oscillation term is of order , so if , then this term converges with a higher rate than and we may assume that it is negligible.

### 3.2. Derivation of the error estimator

Before we derive the error estimator, we first write in terms of element and face distributions. For every , we can write

 ⟨jh,w⟩ =⟨∇×Hh,w⟩ =(Hh,∇×w)Ω =∑T∈Th[(∇×Hh,w)T+(Hh,^nT×w)∂T] =∑T∈Th[(∇×Hh,w)T+(Hh,^nT×w)∂T∖∂Ω] =∑T∈Th[(∇×Hh,w)T+(−^nT×Hh,w)∂T∖∂Ω] =∑T∈Th(∇×Hh,w)T+∑f∈FIh(−[[Hh]]t,w)f =:∑T∈Th(jh,T,w)T+∑f∈FIh(jh,f,w)f,

where denotes the application of a distribution to a function, denotes the set of all internal faces, denotes the normal unit vector to pointing outward of , and denotes the tangential jump operator, with , , and and the two adjacent elements of .

Since and is piecewise constant, we have that . Therefore, if and if , and , where is a normal unit vector of and is given by

 Dk(f):= {u∈L2(f)|u(x)=^nf×(v(x)+xw(x)) for some v∈Pk−1(f)3,w∈Pk−1(f)}.

In other words, can be represented by functions on the elements and face distributions on the internal faces.

We define and can write

 (9) ⟨jΔ,w⟩ =∑T∈Th(jΔT,w)T+∑f∈FIh(jΔf,w)f ∀w∈C∞0(Ω)3∪Rk,0(Th),

where

 (10) jΔT:=j|T−jh,T=j|T−∇×Hh|TandjΔf:=−jh,f=[[Hh]]t|f.

We look for a solution of (7) of the form , with and and where denotes the element-wise gradient operator. The term will take care of the element distributions of and the term will take care of the remaining face distributions.

In the following, we firstly describe how to compute in Section 3.2.1 and characterize the remainder in Section 3.2.2. We then describe how to compute the jumps of on internal faces in Section 3.2.3 and explain how to reconstruct from its jumps in Section 3.2.4.

#### 3.2.1. Computation of ^HΔ

We compute by solving the local problems

 (11a) ∇×^HΔ|T =jΔT, (11b) (μ^HΔ,∇ψ)T =0 ∀ψ∈Pk′(T),

for each element . This problem is well-defined and has a unique solution due to the discrete exact sequence property

 Pk′(T)\xlongrightarrow∇Rk′(T)\xlongrightarrow∇× Dk′(T)\xlongrightarrow∇⋅ Pk′−1(T),

and since and . This last property follows from the fact that due to assumption A2 and due to assumption A1.

#### 3.2.2. Representation of the remainder jΔ−∇×^HΔ

Set

 ^\jΔ :=jΔ−∇×^HΔ.

For every , we can write

 ⟨^\jΔ,w⟩ =⟨jΔ−∇×^HΔ,w⟩ =⟨jΔ,w⟩−⟨∇×^HΔ,w⟩ =⟨jΔ,w⟩−(^HΔ,∇×w)Ω ∑T∈Th(jΔT,w)T+∑f∈FIh(jΔf,w)f−∑T∈Th[(^HΔ,^nT×w)∂T+(∇×^HΔ,w)T] =∑T∈Th[(jΔT−∇×^HΔ,w)T−(^HΔ,^nT×w)∂T]+∑f∈FIh(jΔf,w)f ∑T∈Th[(0,w)T−(^HΔ,^nT×w)∂T∖∂Ω]+∑f∈FIh([[Hh]]t,w)f =∑T∈Th(^nT×^HΔ,w)∂T∖∂Ω+∑f∈FIh([[Hh]]t,w)f =∑f∈FIh([[^HΔ]]t,w)f+([[Hh]]t,w)f =∑f∈FIh([[Hh+^HΔ]]t,w)f,

so

 (12) ⟨^\jΔ,w⟩ =∑f∈FIh([[Hh+^HΔ]]t,w)f=:∑f∈FIh(^\jΔf,w)f

for all . This means that can be represented by only face distributions, and since and , we have that and therefore .

#### 3.2.3. Computation of the jumps of ϕ on internal faces

It now remains to find a such that

 ∇×∇hϕ =^\jΔ.

For every , we can write

 ⟨∇×∇hϕ,w⟩ =(∇hϕ,∇×w)Ω =∑T∈Th[(∇×∇ϕ,w)T+(∇ϕ,^nT×w)∂T] =∑T∈Th[(0,w)T+(∇ϕ,^nT×w)∂T∖∂Ω] =∑T∈Th−(^nT×∇ϕ,w)∂T∖∂Ω

Therefore, we need to find a such that

 −[[∇ϕ]]t|f =^\jΔf ∀f∈FIh.

To do this, we define, for each internal face , the scalar jump with , two orthogonal unit tangent vectors and such that , differential operators , and the gradient operator restricted to the face: . We can then write

 −[[∇ϕ]]t|f =−(^n+×∇ϕ++^n−×∇ϕ−)|f =−(^nf×∇ϕ+−^nf×∇ϕ−)|f =−(^nf×∇fϕ+−^nf×∇fϕ−)|f =−^nf×∇f[[ϕ]]f

for all . We therefore introduce an auxiliary variable and solve

 (13a) −^nf×∇fλf =^\jΔf, (13b) (λf,1)f =0,

for each , where (13b) is only added to ensure a unique solution. In the next section, we will show the existence of and how to construct a such that for all . Now, we will prove that problem (13) uniquely defines . We start by showing that (13) corresponds to a 2D curl problem on a face. To see this, note that . If we take the inner product of (13a) with and , we obtain

 (14a) ∂t2λf =^\jΔf,t1, (14b) −∂t1λf =^\jΔf,t2,

where , which is equivalent to a 2D curl problem on . To show that (13) is well-posed, we use the discrete exact sequence in 2D:

 R\xlongrightarrow⊂ Pk′(f)\xlongrightarrowcurlfDk′(f)\xlongrightarrow∇⋅ Pk′−1(f),

where . Since , it suffices to show that . To prove this, we use that, for every ,

 ⟨^\jΔ,∇ψ⟩ =⟨jΔ−∇×^HΔ,∇ψ⟩ =⟨∇×(H−Hh−^HΔ),∇ψ⟩ =(H−Hh−^HΔ,∇×∇ψ)Ω =(H−Hh−^HΔ,0)Ω =0.

Then, for every , we can write

 0 =⟨^\jΔ,∇ψ⟩ =∑f∈FIh(^\jΔf,∇fψ)f =∑f∈FIh[(^n∂f⋅^\jΔf,ψ)∂f−(∇f⋅^\jΔf,ψ)f] =∑f∈FIh[(^n∂f⋅^\jΔf,ψ)∂f∖∂Ω−(∇f⋅^\jΔf,ψ)f] =∑f∈FIh∑e:e⊂∂f∖∂Ω(^ne,f⋅^\jΔf,ψ)e−∑f∈FIh(∇f⋅^\jΔf,ψ)f =∑e∈EIh∑f:∂f⊃e(^ne,f⋅^\jΔf,ψ)e−∑f∈FIh(∇f⋅^\jΔf,ψ)f =∑e∈EIh⎛⎝