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.

  • 2
    $\begingroup$ You need to distinguish between the residual of the linear system, $\mathbf b - A \mathbf x$ and the residual of the PDE. You need the latter, and you've already stated how it looks: it is a function of the spatial variable $x$, and it can be evaluated at quadrature points. $\endgroup$ Jun 17, 2021 at 15:07
  • $\begingroup$ @WolfgangBangerth so does that mean that using the solution uh, we compute Δ uh (which is zero, if we assume linear elements, as Abdullah Ali Sivas said in his reply) and hence the residual ends up being just the jump terms across faces + the source? $\endgroup$
    – gnikit
    Jun 27, 2021 at 17:14
  • 1
    $\begingroup$ Yes, exactly. On triangles with linear elements, the cell residual is zero and only the jump residual is nonzero. $\endgroup$ Jun 29, 2021 at 3:39

1 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$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.