Residual-type a posteriori error estimator

[VexGH]

Hello,

Assuming that I have the 2D diffusion-reaction problem with pure Neumann b.c., is there a way to obtain a residual type estimator for H1 norm, in order to use it in an adaptive algorithm? Namely, I would like to compute the following quantity for every element K on the mesh,

$$ \eta_{1, K}=\left\{h_K^2\left|f+\Delta u_h-u_h\right|_K^2+\frac{1}{2} \sum_{S \in \partial K \cap S_{int}} h_S\left|\left[ \nabla u_h\cdot \vec{n}\right]\right|_S^2+\sum_{S\in \partial K\cap \partial \Omega} h_S\left|g- \nabla u_h\cdot \vec{n}\right|_S^2\right\}^{\frac{1}{2}}. $$

Here $f$ and $g$ are the right hand side on the domain and the boundary, respectively, while $S$ is an edge and $S_{int}$ is the set of all interior edges.

Best regards,
Vex

[ksagiyam]

I think it is possible. If I understand your problem correctly, you'd like to compute that value on each cell and store that value in a DG0 field?
The cell integral on the r.h.s. looks straightforward. The facet integrals might need to be first projected onto a HDivTrace0 space, and then added to the DG0 field as needed by, e.g., using op2.ParLoop directly.

Other people might know a better way, but, it would help if you could make a small toy problem illustrating what you'd like to do.

[VexGH]

Hello,

Thank you for your answer. I would like to make an adaptive algorithm based on the above residual a posteriori error estimator
and for that reason I need to compute these norms for each cell on the mesh. The DG0 field was indeed my initial thought, but I
did not project the facet integrals onto a HDivTrace0 space. I also did not understand the op2.ParLoop part (I am new to firedrake).

I attach a toy problem with a code, namely, the 2D Poisson equation with pure Robin boundary conditions in order to avoid the
compatibility condition.

The code (GitRobinEst.txt) contains the solution and the H1-error for the above the problem. I made some attemts to compute the
corresponding global and local a posteriori error estimator but I am not sure that any of them is correct. You will find some of those
attempts at the end of the file.

Best Regards,
Vex
GitRobinEst.txt

[ksagiyam]

The approach you took for the 2nd attempt looks good to me and better than what I initially had in my mind.

If you want to take sqrt before sum as in the definition, you will need to modify your code as in the following:

#2nd Alt DG+Assemble (global+local)
DG0 = FunctionSpace(mesh, "DG", 0)
w = TestFunction(DG0)
residual = h**2*w*(div(grad(uh))+fd)**2*dx+avg(w)*avg(h)*jump(grad(uh),n)**2*dS+w*h*(gb-dot(grad(uh),n)-rc*uh)**2*ds
eejv=assemble(residual)
eta1_K = Function(DG0, val=numpy.sqrt(eejv.dat.data_ro))
your_global_norm = numpy.sum(eta1_K.dat.data_ro)
print(your_global_norm)

where eta1_K is a DG0 function that contains the $\eta_{1,K}$ value on each cell.