Computing the residual in a Dual Weighted Residual (DWR) method

Question

I am in the process if computing the Dual-Weighted Residual (DWR) for a linear PDE with a linear functional but I am struggling with the residual part of the calculation.

For example suppose we want to formulate the DWR for a Poisson's equation then the error for a linear functional $J(\cdot)$ is given by

$$ J(e) = \sum_{K} \{ (R_h, z-\varphi_h)_K + (r_h,z-\varphi_h)_{\partial K} \}, $$

and the upper error bound is given by

$$ \eta_K = |(R_h, z-\varphi_h)_K + (r_h,z-\varphi_h)_{\partial K}| . $$

Finally the a posteriori error estimate can be written as

\begin{equation} J(e) \leq \eta = \sum_{K} \eta_{K}. \end{equation}

For the Poisson's equation the element local residuals are

\begin{align} R(u_h)_K & = f + \Delta u_h, \\ r(u_k)_{\Gamma} & = \begin{cases} 1/2 n\cdot[\nabla u_h] && , \text{ if } \Gamma \subset {\partial K} / {\partial \Omega}\\ 0 &&, \text{ if } \Gamma \subset \partial \Omega. \end{cases} \end{align}

for a more detailed explanation see Equations (3.17) and (3.18) from

R. Becker and R. Rannacher, “An optimal control approach to a posteriori error estimation in finite element methods,” Acta Numer., vol. 10, no. 2001, pp. 1–102, 2001.

---

My question is how would one go about computing $R(u_h)_K$ (and $r(u_h)_K$)?

We do have the discretised linear system of the weak form for the Laplacian operator, which is how we obtained $u_h$ in the first place. Wouldn't doing $\mathbf{R} = b - \mathbf{A}x$ for element $K$ be zero since this is a Weighted Residual method or am I missing/mixing something?

NOTE: I am aware of the Galerkin orthogonality property and that something needs to be done when the dual problem's solution error is approximated otherwise the DWR will be zero. Chapter 4 from

W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, 1st ed. Birkhäuser Basel, 2003.`

provides some excellent ideas on the matter.

Answer 1

Take my advice with a bit of salt, as I am not an expert on adaptive FEM. I don't have access to the papers, so I am not sure if the following is how they do it, but it is how I would implement it.

$u_h$ is a polynomial inside each element, so you can take the second derivative of it easily. If the basis consists only up to linear polynomials, $R(u_h)_K$ will be simply equal to $f$ with this approach. I don't know if it is desirable to you.

Similarly, the jump $[\nabla u_h]$ over a face $F$ can be computed as $\nabla u_h^+ - \nabla u_h^-$ where $u_h^+$ and $u_h^-$ are restriction of $u_h$ to the elements $+$ and $-$ which share the face $F$.