Ctrl+K

Equilibrated Residual Error Estimates

29. Equilibrated Residual Error Estimates#

29.1. General framework#

Equilibrated residual error estimators provide upper bounds for the discretization error in energy norm without any generic constant. We consider the standard problem: find \(u \in V := H_0^1(\Omega)\) such that

\[ \int_\Omega \lambda \nabla u \cdot \nabla v = \int_\Omega f v \qquad \forall \, v \in V \]

The left hand side defines the bilinear-form \(A(\cdot, \cdot)\), the right hand side the linear-form \(f(\cdot)\). We define a finite element sub-space \(V_h \subset V\) of order \(k\), and the finite element solution

\[ \text{find } u_h \in V_h: \quad A(u_h, v_h) = f(v_h) \qquad \forall \, v_h \in V_h. \]

We assume that \(f\) is element-wise polynomial of order \(k-1\), and \(\lambda\) is element-wise constant and positive.

The residual \(r(\cdot) \in V^*\) is

\[ r(v) = f(v) - A(u_h, v) \qquad v \in V \]

Since

\[ \| u - u_h \|_A = \sup_{v \in V} \frac{A(u-u_h, v)}{\| v \|_A} = \sup_{v \in V} \frac { r(v) }{\| v \|_A}, \]

we aim in estimating \(\|r \|\) in the norm dual to \(\| \cdot \|_A\), which is essentially the \(H^{-1}\)-norm. In general, the direct evaluation of this norm is not feasible. Using the structure of the problem, we can represent the residual as

\[ r(v) = \sum_{T \in \mathcal{T}} \int_T r_T v + \sum_{E \in \mathcal{E}} \int_E r_E v, \]

where \(r_T\) and \(r_E\) are given given as

\[ r_T = f_T + \operatorname{div} \lambda_T \nabla u_{h|T} \qquad \text{and} \qquad r_E = \left[ \lambda \frac{\partial u_h} {\partial n} \right]_E \]

The element-residual \(r_T\) is a polynomial of order \(k-1\) on the element \(T\), and the edge residual (the normal jump) is a polynomial of order \(k-1\) on the edge \(E\).

The residual error estimator estimates the residual in terms of weighted \(L_2\)-norms:

\[ \| r \|^2 \simeq \eta^{res}(u_h, f)^2 := \sum_T \frac{ h_T^2}{\lambda_T} \| r_T \|_{L_2(T)}^2 + \sum_E \frac{ h_E }{\lambda_E} \| r_E \|_{L_2(E)}^2 \]

Here, \(\lambda_E\) is some averaging of the coefficients on the two elements containing the edge \(E\). The equivalence holds with constants depending on the shape of elements, the relative jump of the coefficient, and the polynomial order \(k\).

The equilibrated residual error estimator \(\eta^{er}\) is defined in terms of the same data \(r_T\) and \(r_E\). It satisfies

\[\begin{align*} \| u - u_h \|_A & \leq \eta^{er} \qquad \text{reliable with constant 1} \\ \| u - u_h \|_A & \geq c \, \eta^{er} \qquad \text{efficient with a generic constant $c$} \\ \end{align*}\]

The lower bound depends on the shape of elements and the coefficient \(\lambda\), but is robust with respect to the polynomial order \(k\).

The main idea is the following: Instead of calculating the \(H^{-1}\)-norm of \(r\), we compute a lifting \(\sigma^\Delta\) such that \(\operatorname{div} \sigma^\Delta = r\), and calculate the \(L_2\)-norm of \(\sigma^\Delta\). Since \(r\) is not a regular function, the equation must be posed in distributional form:

\[ \int_\Omega \sigma^\Delta \cdot \nabla \varphi = - r(\varphi) \qquad \forall \, \varphi \in V \]

Then, the residual can be estimated without envolving any generic constant:

\[\begin{eqnarray*} \| r \|_{A^\ast} & = & \sup_{v \in V} \frac{r(v)}{\| v \|_A} = \sup_v \frac{\int \sigma^\Delta \cdot \nabla v }{\|v \|_A} \\ &= & \sup_v \frac{\int \lambda^{-1/2} \sigma^\Delta \cdot \lambda^{1/2} \nabla v }{\|v \|_A} \leq \sup_v \frac{\sqrt{\int \lambda^{-1} |\sigma^\Delta|^2} \sqrt{ \int \lambda |\nabla v|^2 } }{\|v \|_A} = \| \sigma^\Delta \|_{L_2, 1/\lambda} \end{eqnarray*}\]

The norm \(\| \sigma^\Delta \| := \int \lambda^{-1} | \sigma^\Delta|^2\) can be evaluated easily.

Remark: The flux-postprocessing \(\sigma := \lambda \nabla u_h + \sigma^\Delta\) provides a flux \(\sigma \in H(\operatorname{div})\) such that \(\operatorname{div} \sigma = f\), i.e. the flux is in exact equilibrium with the source \(f\). Thus the name.

29.2. Computation of the lifting \(\| \sigma^\Delta \|\)#

The residual is a functional of the form

\[ r(v) = \sum_T (r_T, v)_{L_2(T)} + \sum_E (r_E, v)_{L_2(E)}, \]

where \(r_T\) and \(r_E\) are polynomials of order \(k-1\). We search for \(\sigma^\Delta\) which is element-wise a vector-valued polynomial of order \(k\), and not continuous across edges. Element-wise integration by parts gives

\[ \int_\Omega \sigma \cdot \nabla \varphi = -\sum_T \int_T \operatorname{div} \sigma_{|T} \varphi + \sum_E \int_E [\sigma \cdot n]_E \varphi. \]

Thus \(\operatorname{div} \sigma = r\) in distributional sense reads as

\[ \operatorname{div} \sigma_{|T} = r_T \qquad \text{and} \qquad [\sigma \cdot n]_E = -r_E \]

for all elements \(T\) and edges \(E\). We could now pose the problem

\[ \min_{ \sigma \in P^k({\mathcal T})^2 \atop \operatorname{div} \sigma = r} \| \sigma \|_{L_2, 1/\lambda} \]

We minimize the weighted-\(L_2\) norm since we want to find the smallest possible upper bound for the error. This is already a computable approach. But, the problem is global, and its solution is of comparable cost as the solution of the original finite element system. The existence of a \(\sigma\) such that \(\operatorname{div} \, \sigma = r\) also needs a proof.

We want to localize the construction of the flux. Local problems are associated with vertex-patches \(\omega_V = \cup_{T : V \in T} T\). We proceed in two steps:

  • localization of the residual: \(r = \sum_V r^V\)

  • local liftings: find \(\sigma^V\) such that \(\operatorname{div} \sigma^V = r^V\) on the vertex patch. Then, for \(\sigma := \sum \sigma^V\) there holds \(\operatorname{div} \sigma = r\)

The localization is given by multiplication of the \(P^1\) vertex basis functions (hat-functions) \(\phi_V\):

\[ r^V(v) := r( \phi_V v) \]

Since \(\sum_V \phi_V = 1\), there holds \(\sum r^V(\cdot) = r(\cdot)\). The localized residual has the same structure of element and edge terms:

\[ r^V(v) = \sum_{T \subset \omega_V} (r_T^V, v)_{L_2(T)} + \sum_{E \subset \omega_V} (r_E^V, v)_{L_2(E)}, \]

with

\[ r_T^V = \phi_V r_T \qquad \text{and} \qquad r_E^V = \phi_V r_E \]

The local residual vanishes on constants on the patch:

\[ r^V(1) = r(\phi_V \, 1) = A(u-u_h, \phi_V) = 0 \]

The last equality follows from the Galerkin-orthogonality.

We give an explicit construction of the lifting \(\sigma^V\) in terms of the Brezzi-Douglas-Marini (BDM) element. The \(k^{th}\) order BDM element on a triangle is given by \(V_T = [P^k]^2\) and the degrees of freedom:

  • \(\int_E \sigma\cdot n \, q_i\) with \(q_i\) a basis for \( P^k(E)\)

  • \(\int_T \operatorname{div} \sigma \, q_i\) with \(q_i\) a basis for \( P^{k-1}(T) \cap L_2^0(T)\)

  • \(\int_T \sigma \cdot \operatorname{curl} q_i\) with \(q_i\) a basis for \( P_0^{k+1}(T)\)

Exercise: Show that these dofs are unisolvent. Count dimensions, and prove that \([ \forall i: \psi_i(\sigma) =0 ] \Rightarrow \sigma = 0\).

Now, we give an explicit construction of equilibrated fluxes on a vertex patch. Label elements \(T_1, T_2, \ldots T_n\) in a counter-clock-wise order. Edge \(E_i\) is the common edge between triangle \(T_{i-1}\) and \(T_i\) (with identifying \(T_0 = T_n\)). We define \(\sigma\) by specifying the dofs of the BDM element:

  • Start on \(T_1\). We set \(\sigma_n = -r^V_{E_1}\) on edge \(E_1\). On the edge on the patch-boundary we set \(\sigma_n = 0\), and on \(E_2\) we set \(\sigma_n = const\) such that \(\int_{\partial T_1} \sigma_n = \int_{T_1} r^V_T\). We use the dofs of type (ii) to specify \(\int_T \operatorname{div} \sigma \, q = \int_T r_T^V q \; \forall q \in P^{k-1} \cap L_2^0(T)\). Together with get \(\operatorname{div} \sigma = r_T\). Dofs of type (iii) are not needed, and set 0. There holds $\( \int_{E_2} \sigma_n = \int_{T_1} r^V_T - \int_{E_1} \sigma_n = \int_{T_1} r^V_{T_1} + \int_{E_1} r^V_{E_1} \)$

  • Continue with element \(T_2\). On edge \(E_2\) common with \(T_1\) set \(\sigma_n\) such that \([\sigma \cdot n]_{E_2} = r_{E_2}\). Otherwise, proceed as on \(T_1\). Thus $\( \int_{E_3} \sigma_n = \int_{T_1} r^V_{T_1} + \int_{E_1} r^V_{E_1} +\int_{T_2} r^V_{T_2} + \int_{E_2} r^V_{E_2} \)$

  • Continue to element \(T_n\). Observe that on \(T_n\): $\( \int_{E_1} \sigma_n = \sum_{i=1}^n \int_{T_i} r^V_{T_i} + \sum_{i=1}^n \int_{E_i} r^V_{E_i} = 0, \)\( which follows from \)r^V(1) = 0\(. Thus, also \)[\sigma \cdot n]{E_1} = r{E_1}^V$ is satisfied.

This explicit construction proves the existence of an equilibrated flux. Instead of this explicit construction, one may solve a local constrained optimization problem

\[ \min_{\sigma^V : \operatorname{div} \sigma^V = r^V} \| \sigma \|_{L_2, \lambda^{-1}} \]

This applies also for 3D. Furthoer notes

  • mixed boundary conditions are possible

  • the efficiency for the h-FEM is shown by scaling arguments, and equivalence to the residual error estimator

  • efficiency is also proven to be robust with respect to polynomial order \(k\), examples show overestimation less than 1.5

Literature:

  • D. Braess and J. Schöberl: Equilibrated Residual Error Estimator for Maxwell’s Equations. Mathematics of Computation, Vol 77(262), 651-672, 2008

  • D. Braess, V. Pillwein and J. Schöberl: Equilibrated Residual Error Estimates are p-Robust. Computer Methods in Applied Mechanics and Engineering. Vol 198, 1189-1197, 2009

29.3. Equilibration in NGSolve#

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *

r1 = Rectangle(1,1).Face()
r1.edges.Max(X).name='right'
r1.edges.Min(X).name='left'
r1.edges.Max(Y).name='top'
r1.edges.Min(Y).name='bottom'
r1.faces.name='air'
r2 = MoveTo(0.6,0.6).Rectangle(0.2,0.2).Face()
r2.faces.name='source'

r1 -= r2
shape = Glue( [r1,r2] )

mesh = Mesh(OCCGeometry(shape,dim=2).GenerateMesh(maxh=0.1))
Draw (mesh);

order = 3
order_equ = 3
dirichlet="left|bottom"
fes = H1(mesh, order=order, dirichlet=dirichlet)
u,v = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
source = mesh.MaterialCF( { "source" : x*y }, default=0 )
f = LinearForm(source*v*dx).Assemble()

gfu = GridFunction(fes, name="solution")
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu);

construct \(\sigma\) such that \(\operatorname{div} \sigma = f + \Delta u_h\) from local contributions \(\sigma = \sum \sigma^V\):

\[ \min_{\operatorname{div} \sigma^V = \phi_V (f + \Delta u_h)} \| \sigma^V \|_0^2 \]

where \(\operatorname{div} \sigma^V = \phi_V (f + \Delta u_h)\) is understood in distributional sense, what is

\[\begin{eqnarray*} \operatorname{div}_T \sigma^V & = & \phi_V (f + \Delta_T u_h) \\ [\sigma_n^V] & = & \phi_V[\partial_n u_h] \end{eqnarray*}\]
# variational formulation for patch correction problems
Xsigma = Discontinuous (HDiv(mesh, order=order_equ))
Xw     = L2(mesh,order=order_equ-1)
Xwf    = FacetFESpace(mesh, order=order_equ)
Xlam   = NumberSpace(mesh)
X = Xsigma*Xw*Xwf*Xlam 

sigma, w, wf, lam = X.TrialFunction()
tau,   v, vf, mu  = X.TestFunction()

n = specialcf.normal(mesh.dim)
bfequ = (tau*sigma + w*div(tau) + v*div(sigma) + w*mu + lam*v) * dx 
bfequ += (-sigma*n*vf-tau*n*wf)*dx(element_boundary=True)  
bfequ += (-1e10*wf.Trace()*vf.Trace())*ds(definedon=mesh.Boundaries(dirichlet)) # penalty for Dirichlet
bfequ += (-1e10*mu*lam)*ds(definedon=mesh.Boundaries(dirichlet))  # no constraint at Dir bnd

# the residual: element-term + edge term:
lfequ = (source+Trace(gfu.Operator("hesse")))*v*dx 
lfequ += (-grad(gfu)*n*vf)*dx(element_boundary=True)

gfequ = GridFunction(X)

with TaskManager():
    PatchwiseSolve (bf=bfequ, lf=lfequ, gf=gfequ)

fesflux = HDiv(mesh, order=order_equ)
equflux = GridFunction(fesflux)
equflux.Set (grad(gfu)-gfequ.components[0])

Draw (equflux, mesh);

Check equilibrium:

Draw (div(equflux)+source, mesh, "residuum");
Contents