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");

Equilibrated Residual Error Estimates

Contents

29. Equilibrated Residual Error Estimates#

29.1. General framework#

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

29.3. Equilibration in NGSolve#