A posteriori error estimation based on the equilibrated flux recontruction

The Laplace problem

(In the following, only some extractions of the full script will be shown to help understand the implementation work. You can download the full script at the end of this page.)

Let $\Omega\subset \mathbb{R}^2$ be an open bounded polygonal domain with the boundary donoted by $\partial \Omega$. Let us consider here the Laplace equation: for $f\in L^2(\Omega)$, find $u: \Omega \rightarrow \mathbb{R}$ such that \begin{align} -\Delta u = f &\qquad \mbox{ in }\Omega, \\ u = 0 & \qquad \mbox{ on } \partial\Omega. \end{align} The variational formulation is defined by: Find $u\in H^1_0(\Omega)$ such that $$(\nabla u, \nabla v) = (f,v) \qquad \forall v\in H^1_0(\Omega).$$ We denote by $\mathcal{T}_h$ the mesh used in the following. $\mathcal{T}_h$ is a simplicial partition of the polytope $\Omega$, i.e. $\cup_{K\in\mathcal{T}_h} \overline{K} = \overline{\Omega}$, any $K\in\mathcal{T}_h$ is a closed simplex (a triangle in 2D), and the intersection of two different simplices is either empty (a vertex, an edge in 2D).

Here we consider the computational domain as a square, with an uniform mesh. The number of the node on each edge is denoted by numnode.
// define the mesh
int numnode = 80;
mesh Th = square(numnode,numnode);

For this example, we consider the $\mathbb{P}_1$-conforming finite element method, that means we shall be using $V_h\subset H^1_0(\Omega)$ $$V_h=\mathbb{P}_1(\mathcal{T}_h)\cap H^1_0(\Omega) = \left\{v_h\in H^1_0(\Omega); v_h\mid_K\in \mathbb{P}_1(K), \qquad K\in \mathcal{T}_h\right\}.$$
// functional space for solution
fespace Vh(Th,P1);

Using the finite element method, the discrete variational formulation is given by: Find $u_h\in V_h$ such that $$(\nabla u_h, \nabla v_h) = (f,v_h) \qquad \forall v\in V_h.$$
// boundary condition
func gd=x*(x-1)*y*(y-1);

func f = -2*(x*x+y*y)+2*(x+y);

// variational formulation

// solve the laplace problem exactly
matrix A=a(Vh,Vh,solver=CG);
real[int] b=a(0,Vh);
uh[]=A^-1*b;

We call the equilibrated flux reconstruction any function $\sigma_h$ constructed from $u_h$ which satisfies $\sigma_h\in \mathbf{H}(\rm div,\Omega)$, $$(\nabla\cdot \sigma_h, 1)_K= (f, 1)_K, \qquad \forall K\in \mathcal{T}_h$$ For any $K\in \mathcal{T}_h$, define the residual estimator by: $$\eta_{{\rm R},K} := \frac{h_K}{\pi} \|f-\nabla \cdot \sigma_h\|_K,$$ the flux estimator by $$\eta_{{\rm F},K} := \|\nabla u_h+ \sigma_h\|_K.$$ Then we have $$\|\nabla (u-u_h)\|^2 \leq \sum_{K\in \mathcal{T}_h} (\eta_{{\rm F},K}+\eta_{{\rm R},K} )^2$$ To construct the equilibrated flux, we should solve the local Neumann mixed finite element problem, some more technical details can be found here (page 54--56).

To obtain these equilibrated flux, we need to solve the local mixed finite element problem by the following four steps.

1) we construct the "patches" associated to each node in the mesh.

mesh ThaI=Th;
mesh[int] Tha(Th.nv);
matrix[int] Aa(Th.nv);

fespace Pa(ThaI,P1dc);
for(int i=0; i < Th.nv; ++i)
{
phia[][i]=1; // phia
Tha[i] = trunc (Th, phia > 0,label=10);
if (Th(i).label==0) Tha[i]=change(Tha[i],flabel=10);
phia[][i]=0; // phia
}



2) For the local problem, we need to construct the matrix

fespace VPa(ThaI,[RT1,P1dc]);

real epsregua = numnode*numnode*1e-10;   //regularization for local problems

for( int i=0; i < Th.nv ;++i)
{
//  to bluild the basic function
phia[][i]=1; // phia
ThaI = Tha[i];
VPa [s1,s2,ra],[v1,v2,q];
//matrix
varf a([s1,s2,ra],[v1,v2,q]) = int2d(ThaI)( [s1,s2]'*[v1,v2] - ra*Div(v1,v2) - Div(s1,s2)*q - epsregua*ra*q )+ on(10, s1=0,s2=0);
matrix Art=a(VPa,VPa);
Aa[i]= Art;
set(Aa[i],solver=sparsesolver);
phia[][i]=0;
}

Here we add one regularization term as $-\int_{ThaI} epsregua*ra*q$ to ensure the average value of $ra$ is 0.

3) We need also the right-hand-side of the local problems, and we solve them.
//for loop computing right hand side and solve the local problems
fespace Pa(ThaI,P1dc);
fespace VRT(Th,RT1);

VRT [sigma1,sigma2]=[0,0];
for(int i=0; i < Th.nv ;++i)
{
phia[][i]=1; // phia
ThaI = Tha[i];
VPa [s1,s2,ra],[v1,v2,q];
//right hand side
+ on(10,uv1 = 0, uv2 = 0);
VPa [F,F1,F2] ; F[] = l(0,VPa);
Fa[i].resize(F[].n);
Fa[i]=F[];
s1[] = Aa[i]^-1*Fa[i];
real[int] so(s1[].n);
so=sigma1[](I[i]);
s1[]= s1[].*epsI[i];
s1[] += so;
sigma1[](I[i]) = s1[];
phia[][i]=0;
}


4) with this equilibrated flux, we can compute our estimator
// computation for the estimator

+ int2d(Th)(chiK*(hTriangle*hTriangle/pi/pi* (f-Div(sigma1,sigma2))*(f-Div(sigma1,sigma2)) ));


In order to showcase the performance of the estimator, we consider an example with an analytic solution,

// case with an analytic solution, and its derivatives [dx(u),dy(u)]
func uExact = x*(x-1)*y*(y-1);
func dxuExact=(2*x-1)*y*(y-1);
func dyuExact=(2*y-1)*x*(x-1);
// boundary condition
func gd=x*(x-1)*y*(y-1);
func f = -2*(x*x+y*y)+2*(x+y);


We compute the exact error (discretization error)
// computation for the exact error and map of error distribution
fespace Wh(Th,P0); // constant in each element
Wh mapErrorTotal;
mapErrorTotal[]=indicatorErrorTotal(0,Wh);
real GlobalError=sqrt(mapErrorTotal[].sum);
mapErrorTotal=sqrt(mapErrorTotal);
plot(mapErrorTotal, ps="MapDiscretisationError.eps",value=1, fill=1,cmm="DiscretisationError");

We show here the analytic solution ($u$) and the solution obtained by the finite element method ($u_h$)

 Analytic solution $u$ Numerical solution $u_h$

as well as the distribution of the discretization error $\|u-u_h\|_K$ and the estimator $\eta_{{\rm F},K}+\eta_{{\rm R},K}$

 Discretization error Discretization estimator

With five levels of uniform mesh refinement, the effectivity index (the ratio between the total estimator and the total error) is given as follows:
 Effectivity index