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.