An a posteriori error bound for discontinuous Galerkin approximations of convection–diffusion problems

An a posteriori error bound for discontinuous Galerkin approximations of convection–diffusion... Abstract An a posteriori bound for the error measured in the discontinuous energy norm for a discontinuous Galerkin (dG) discretization of a linear one-dimensional stationary convection–diffusion–reaction problem with essential boundary conditions is presented. The proof is based on a conforming recovery operator inspired from a posteriori error bounds for the dG method for first-order hyperbolic problems. As such, the bound remains valid in the singular limit of vanishing diffusion. Detailed numerical experiments demonstrate the independence of the quality of the a posteriori bound with respect to the Péclet number in the standard dG-energy norm, as well as with respect to the viscosity parameter. 1. Introduction The interaction between convection and diffusion in many physical systems often gives rise to multiscale solution behaviour, when the convection is the dominant phenomenon. These are typically manifested as boundary and/or interior layers or even as discontinuities, in cases of diffusion-degenerate problems. The stable discretization for convection–diffusion problems in the contexts of finite difference and finite element methods is now relatively well understood. However, such lower-dimensional solution features require carefully selected variable mesh resolution across the computational domain to be efficiently resolved. The location of such multiscale features may not be a priori available. This motivates the need for adaptive algorithms. In the context of finite element methods, there exists a growing literature on the derivation of adaptive algorithms for (steady-state and transient) convection–diffusion problems, based on local error estimation, either in an ad hoc fashion or via reliable a posteriori bounds. A key question in the respective literature is the independence (robustness in the terminology of most references) of the derived a posteriori estimators with respect to the, so-called, Péclet number, i.e., the ratio between the convection and the diffusion magnitudes. To achieve such error control with upper and lower a posteriori bounds, several approaches have been considered, with the most popular being the augmentation of the energy norm by a negative/fractional Sobolev norm of the skew-symmetric part of the differential operator (Verfürth, 2005; Sangalli, 2008; Schötzau & Zhu, 2009; Ern et al., 2010; Creusé & Nicaise, 2013); we refer to Sangalli (2008) for an insightful exposition of the differential operator’s stability properties with respect to the Péclet number. We also mention the use of subgrid problems (Sangalli, 2001), the recent adaptive variational stabilization methods (Cohen et al., 2012; Dahmen et al., 2012; Ainsworth et al., 2013), focusing on fully computable a posteriori bounds for the standard energy norm. Generally speaking, the aforementioned works (with the exception of Cohen et al., 2012; Dahmen et al., 2012, which introduce new stabilized numerical methods aiming at sharp error control) are based on adapting the a posteriori error analysis for the purely elliptic problem to the convection–diffusion(–reaction) problem yielding error estimates controlling the error with upper and lower bounds. As a result, although the a posteriori error bounds in Verfürth (2005), Sangalli (2008), Ern et al. (2010), Schötzau & Zhu (2009) and Creusé & Nicaise (2013) hold for each fixed positive diffusion parameter, the derived estimators are not applicable when the diffusion vanishes; this limits the possible applicability of such bounds to realistic nonlinear shock computations. Given the very successful computational performance of standard stabilized methods, such as discontinuous Galerkin (dG) methods with various fluxes, we would like to seek an alternative approach to their error control. Our work is motivated by the following issues: (1) whether it is possible to prove a posteriori bounds for singularly perturbed problems of this kind starting from a posteriori error control of the hyperbolic limit problem and (2) whether such a posteriori bounds for errors measured in standard energy-type norms (i.e., norms with respect to which the method’s bilinear forms are coercive) would hold with error bound constants that remain uniformly bounded with respect to vanishing viscosity and the Péclet number. Our answer to both questions is that, under certain conditions, such error bounds are possible. We stress that, in the present context of convection–diffusion problems, we address these questions with regard to the proof of upper bounds and their respective validity at the zero viscosity limit. In addition, we address, to a certain extent, similar questions related to lower bounds and we thoroughly investigate computationally the behaviour of the estimators. We derive a new a posteriori error bound for dG discretizations of a linear stationary convection–diffusion–reaction problem, utilizing ideas from error control for first-order hyperbolic problems. In fact, we show that the answers to the questions (1) and (2) above are positive, provided we have at our disposal a technique that leads to error control of the limiting hyperbolic problem. This can be done to the full generality considering the one-dimensional problem by utilizing the reconstructions in Makridakis & Nochetto (2006). The multidimensional hyperbolic problem was partially addressed in Georgoulis et al. (2014), and its full generality is an important open problem. Thus, we have chosen to present our results in full detail for the one-dimensional problem. It is clear from the analysis that similar results can be obtained in higher dimensions using the ideas in Georgoulis et al. (2014). Since such results are still not available in their full generality, we have chosen to consider numerical experiments to highlight the applicability of the proposed approach in higher dimensions. To this end, we derive a posteriori upper bounds for the dG method for the convection–diffusion–reaction problem, which remain valid in the hyperbolic limit of vanishing diffusion and at the same time remain bounded uniformly with respect to the Péclet number (see Theorems 4.2 and 4.3 and Section 4.2). The upper bounds are complemented with standard lower bounds of the flux jump term whose constant, as expected, degenerates with the Péclet number tending to infinity. To further theoretically justify the computationally observed robustness of the proposed a posteriori estimators with respect to the Péclet number, additional nonstandard lower bounds are also derived using classical asymptotic analysis results (cf. Theorem 5.2). The asymptotically optimal behaviour of the a posteriori error estimator, as well as its superiority on controlling the natural dG-energy norm of the problem, is observed in practice on a number of one-dimensional numerical examples. Moreover, we observe numerically that the quality of the derived a posteriori bounds with respect to the energy norm does not deteriorate with increasing Péclet number (when the viscosity tends to zero), as opposed to the respective estimators from Verfürth (2005), Sangalli (2008), Ern et al. (2010), Schötzau & Zhu (2009) and Creusé & Nicaise (2013). To this end, numerical comparisons with known estimators from the literature are performed. Finally, we consider some two-dimensional numerical experiments to highlight the performance of the estimators in higher dimensions. The remainder of this work is structured as follows. In Section 2, we present the model problem and the dG method for the one-dimensional case. In Section 3, we define and discuss the properties of the reconstruction used in the a posteriori analysis in Section 4. Section 4 contains the main results of this work, providing a posteriori bounds for various regimes of equation coefficients, while in Section 5 lower bounds are discussed. Section 6 contains a number of numerical experiments highlighting, in particular, the viscosity independence of the a posteriori error bounds with respect to the Péclet number, along with use of new a posteriori bounds within adaptive algorithms. 2. Preliminaries 2.1 Model problem and notation We shall make use of standard notation for the $$L^2$$-inner product $$(\cdot,\cdot)_{\omega}$$ on an interval $$\omega\subset \mathbb{R}$$, along with the respective $$L^2$$-norm $$\|{\cdot}\|_{\omega}$$; the subscript $$\omega$$ will be suppressed when $$\omega \equiv I:=(\alpha,\beta)\subset\mathbb{R}$$, the computational interval. The standard Hilbertian Sobolev spaces $$H^1(\omega)$$ and $$H^1_0(\omega)$$ will also be used. In $$I$$, we consider the boundary value problem \begin{equation}\label{pde} -\varepsilon u ''+b u'+cu =f\quad\text{in } I,\qquad u(\alpha) = 0 = u(\beta), \end{equation} (2.1) for $$f\in L^2(I)$$, $$\varepsilon > 0$$, $$0\ne b\in \mathbb{R}$$ and $$c:\bar{I}\to\mathbb{R}$$ a continuous function; in weak form this reads, find $$u\in H^1_0(I)$$, such that \begin{equation}\label{pde_weak} \varepsilon (u ',v')+(b u',v)+(cu,v) =(f,v)\quad\text{for all } v\in H^1_0(I). \end{equation} (2.2) We define the inflow node of $$I$$ by $$\partial_- I:=\alpha$$ if $$b>0$$ and $$\partial_- I:=\beta$$ if $$b<0$$; the respective outflow node is defined by $$\partial_+ I:=\{\alpha,\beta\}\backslash \partial_- I$$. The singular limit case $$\varepsilon =0$$ (together with the respective removal of the outflow boundary condition) will also be discussed; the respective hyperbolic problem in weak form reads, find $$u\in H^1_-(I)$$, such that \begin{equation}\label{pde_weak_hyp} (b u',v)+(cu,v) =(\,f,v)\quad\text{for all } v\in H^1_-(I), \end{equation} (2.3) with $$H^1_-(I):=\{v\in H^1(I): v(\partial_- I)=0\}$$. We note that completely analogous a posteriori bounds to the ones given below can be shown for mildly more general coefficients in (2.1) than those assumed above. We shall refrain from doing so in the interest of simplicity of the presentation, settling for providing comments in this direction at various parts of this work. Let $$\mathcal{T}$$ be a subdivision of the domain $$I$$ into elements $$T_i:=(x_{i-1},x_i)$$ for $$i=1,\dots, N$$, $$x_0:=\alpha$$, $$x_{N}:=\beta$$; each element $$T_i$$ has exactly one inflow node $$\partial_-T_i\in \{x_{i-1},x_i\}$$, depending on the sign of $$b$$. Also let $$\partial_+T_i:=\{x_{i-1},x_i\}\backslash \partial_-T_i$$. For notational convenience, we shall also denote by $$\mathcal{N}_-$$ (resp. $$\mathcal{N}_+$$) the set of indices of all interior nodes $$i=1,\dots, N-1$$, together with the index of the inflow node $$\partial_-I$$ (resp. outflow node $$\partial_+I$$). Further, we consider the corresponding (discontinuous) elementwise polynomial space \[ \mathbb{V}_p:=\left\{v\in L^2(I): v|_{(x_{i-1},x_{i})}\in\mathcal{P}_p(x_{i-1},x_{i})\right\}\!, \] where $$\mathcal{P}_p(x_{i-1},x_{i})$$ is the space of polynomials of degree at most $$p\in\mathbb{N}$$. For $$v\in \mathbb{V}_p$$, $$i=0,\dots,N$$, we define the upwind jump $$\lfloor{v}\rfloor(x_i)$$ across the node $$x_i$$ by $$\lfloor{v}\rfloor(x_i):=\lim_{\delta\to 0} ( v(x_i+b\delta) -v(x_i-b\delta) )$$, $$\delta>0$$, adopting the conventions $$\lfloor{v}\rfloor(\alpha) = v(\alpha)$$, $$\lfloor{v}\rfloor(\beta) = v(\beta)$$, i.e., the values taken from outside the domain $$I$$ are set to be equal to the homogeneous boundary conditions. We also define the average and jump across $$x_i$$ by \[ \{ \kern -1.6mm \{{v}\} \kern -1.6mm \}\ (x_i) := (v(x_i^+) +v(x_i^-))/2\quad\text{ and} \quad [ \kern -.7mm [{v}] \kern -.7mm ](x_i) := v(x_i^-) -v(x_i^+), \] respectively, with $$v(x_i^+)$$ and $$v(x_i^-)$$ denoting the traces from the right and from the left, respectively; we also adopt the conventions $$\{ \kern -1.6mm \{{v}\} \kern -1.6mm \}(\alpha)= v(\alpha)$$, $$[ \kern -.7mm [{v}] \kern -.7mm ](\alpha)= -v(\alpha)$$ and $$\{ \kern -1.6mm \{{v}\} \kern -1.6mm \}(\beta)=[ \kern -.7mm [{v}] \kern -.7mm ](\beta) = v(\beta)$$. Using this notation, we can also define the mesh-size function $$h:\bar{I}\to \mathbb{R}_+$$, given by $$h|_{T_i} = h_i:=x_i-x_{i-1}$$, $$i=1,\dots,N$$, $$h(\alpha) = h_1$$, $$h(\beta) = h_N$$, and $$h(x_i) := \{ \kern -1.6mm \{{h}\} \kern -1.6mm \}$$ for $$i=1,\dots,N-1$$. 2.2 The dG method We consider the (interior penalty) dG method reading, find $$u_h\in \mathbb{V}_p$$, such that \begin{equation}\label{dg_weak} \varepsilon B_{\rm d}(u_h,v_h)+ B_{\rm c}(u_h,v_h) = (\,f,v_h) \end{equation} (2.4) for all $$v_h\in \mathbb{V}_p$$, with \[ \begin{aligned} B_{\rm d}(w,v) := & \sum_{i=1}^N (w',v')_{T_i}-\sum_{i=0}^{N} \left(\{ \kern -1.6mm \{{w'}\} \kern -1.6mm \}[ \kern -.7mm [{v}] \kern -.7mm ]+ \{ \kern -1.6mm \{{ v'}\} \kern -1.6mm \}[ \kern -.7mm [{w}] \kern -.7mm ]-\sigma[ \kern -.7mm [{w}] \kern -.7mm ][ \kern -.7mm [{v}] \kern -.7mm ]\right)(x_i)\\ \end{aligned} \] and \[ B_{\rm c}(w,v):= \sum_{i=1}^N\Big( \big(bw'+cw,v\big)_{T_i} + |b|(\lfloor{w}\rfloor v^+)(\partial_-T_i)\Big), \] with $$v_h^+(\partial_-T_i)$$ denoting the value of $$v_h$$ on the inflow node $$\partial_-T_i$$ from within $$T_i$$, and $$\sigma> 0$$ given by $$\sigma(x_i) = C_{\sigma}/h(x_i)$$, $$i=0,\dots,N$$, for some user-defined constant $$C_{\sigma}>1$$. The corresponding dG-energy norm associated to $$B_{\rm d}$$ is defined by \[ | \kern -.25mm \|{w}| \kern -.25mm \| := \left(\sum_{i=1}^N \|{w'}\|_{T_i}^2 +\sum_{i=0}^{N}(\sigma[ \kern -.7mm [{w}] \kern -.7mm ]^2)(x_i)\right)^{\frac{1}{2}}. \] We also define the dG-energy norm for the method by \begin{equation}\label{en_norm} | \kern -.25mm \|{w}| \kern -.25mm \|_q : = \left(\varepsilon| \kern -.25mm \|{w}| \kern -.25mm \|^2 + \|{qw}\|_{}^2 + \frac{|b|}{2} \sum_{i=0}^{N}[ \kern -.7mm [{w}] \kern -.7mm ]^2(x_i)\right)^{\frac{1}{2}}, \end{equation} (2.5) with the subscript indicating the weight in the $$L^2$$-norm component of the dG-energy norm. The bilinear form $$B_{\rm d}$$ can be shown to be coercive in this norm in $$H^1_0\times H^1_0$$, and also in $$\mathbb{V}_p\times \mathbb{V}_p$$, provided the penalty constant $$C_{\sigma}$$ is chosen sufficiently large; for a proof we refer, e.g., to Arnold (1982) and Houston et al. (2002). For the bilinear form $$B_{\rm c}$$, the identity \[ B_{\rm c}(w,w) = \|{\sqrt{c}w}\|_{}^2 + \frac{|b|}{2} \sum_{i=0}^{N}[ \kern -.7mm [{w}] \kern -.7mm ]^2(x_i) \] holds for $$c\ge0$$ (see, e.g., Johnson & Pitkäranta, 1986; Houston et al., 2002). Weaker assumptions on $$c$$ will be discussed in Section 4.3. We note that the theory presented below also remains valid for nonsymmetric and incomplete interior penalty dG methods; this has not been carried through explicitly to minimize the notational overhead. 3. Reconstruction We begin by defining a reconstruction of the approximate solution, which is closely related to the optimal-order time reconstruction for dG time-stepping schemes, presented in Makridakis & Nochetto (2006). This will be a crucial ingredient in the proof of the a posteriori bound below, enabling its validity up to, and including, the singular limit $$\varepsilon =0$$. Definition 3.1 (Reconstruction). For each $$i=1,\dots,N$$, we define the dG reconstruction $$\hat{u}_h\in \mathbb{V}_{p+1}$$ of $$u_h\in \mathbb{V}_p$$ on each $$T_i$$, by the relations \begin{equation}\label{space_rec_def_formula} \big( (b\hat{u}_h)', v_h\big)_{T_i} = \big( (b u_h)', v_h\big)_{T_i} + |b|(\lfloor{u_h}\rfloor v_h^+)(\partial_-T_i) \end{equation} (3.1) for all $$v_h\in \mathbb{V}_p$$, and $$\hat{u}_h(\partial_-T_i) = u_h^-(\partial_-T_i)$$, with $$u_h^-(\partial_-T_i)$$ denoting the value of $$u_h$$ on the inflow node $$\partial_-T_i$$ from outside $$T_i$$; when $$\partial_-T_i=\partial_-I$$, we set $$\hat{u}_h(\partial_-T_i) = 0$$, i.e., the inflow boundary value. The following lemma shows the well-posedness of this reconstruction. The proof is essentially given in Makridakis & Nochetto (2006, Lemma 2.1) in the context of dG time-stepping methods; it is included here for completeness of the presentation. Lemma 3.2. The dG reconstruction $$\hat{u}_h$$ is uniquely defined and we have \begin{equation}\label{space_rec} \big( (b \hat{u}_h)' + cu_h, v_h\big) = B_{\rm c}(u_h,v_h) \end{equation} (3.2) for all $$v_h\in \mathbb{V}^n$$. Moreover, $$\hat{u}_h$$ is continuous in $$I$$. Proof. Identity (3.2) is evident from Definition 3.1. To show the continuity of the reconstruction for all but the outflow element, we integrate (3.1) by parts and we use the property $$\hat{u}_h(\partial_-T_i) = u_h^-(\partial_-T_i)$$, to obtain \begin{equation}\label{reczero} \int_{x_{i-1}}^{x_{i}} b(\hat{u}_h-u_h) v'_h\,\mathrm{d} x =|b|\left(\left(\hat{u}_h-u_h^+\right) v_h^+\right)(\partial_+T_i). \end{equation} (3.3) Now, setting $$v_h$$ to be a constant on $$T_i$$, we deduce $$\hat{u}_h(\partial_+T_i)=u_h^+(\partial_+T_i)$$, with $$u_h^+(\partial_+T_i)$$ denoting the value of $$u_h$$ at the outflow node $$(\partial_+T_i)$$ taken from within $$T_i$$. This implies continuity at all interior nodes $$x_i$$, $$i=1,\dots,N-1$$ and that the left-hand side of (3.3) is equal to zero for all $$v_h'\in\mathcal{P}_{p-1}(x_{i-1},x_i)$$. Hence, $$\pi_{p-1}(\hat{u}_h-u_h)=0$$ on $$(x_{i-1},x_i)$$, with $$\pi_{q}$$ denoting the orthogonal $$L^2$$-projection onto $$\mathbb{V}_q$$. This, together with the exactness at the interval end points, implies the uniqueness of the reconstruction on $$(x_{i-1},x_i)$$. Finally, (3.2) follows by summing for all $$i=1,\dots,N$$. □ A crucial stability property of the above reconstruction is the content of the following result. Lemma 3.3. For each element $$T_i$$, $$i=1,\dots,N$$, we have \[ \left \|{(\hat{u}_h-u_h)^{(s)} }\right\|_{T_i} \le C h_i^{1/2-s}|[ \kern -.7mm [{u_h}] \kern -.7mm ](\partial_-T_i)|, \] $$s=0,1$$, for a $$C>0$$ constant, independent of $$h_i$$ and of $$u_h\in\mathbb{V}_p$$, but dependent on the elemental polynomial degree $$p$$. Proof. The proof for $$s=0$$ follows completely analogously to the one of Makridakis & Nochetto (2006, Lemma 2.2) and is, therefore, omitted. The proof for $$s=1$$ then follows using the standard inverse estimate $$\|{(\hat{u}_h-u_h)' }\|_{T_i} \le Ch_i^{-1} \|{\hat{u}_h-u_h }\|_{T_i}$$ along with the bound for $$s=0$$. □ The dependence of the constant on the polynomial degree $$p$$ in Lemma 3.3 has been exactly understood (Schötzau & Wihler, 2010) but we refrain from retaining this dependence below in the interest of simplicity of presentation. Moreover, it is possible to extend the above to nonconstant $$b$$, provided $$b$$ does not change sign inside an element; this will be considered in detail elsewhere. Lemma 3.4. Let $$\bar{q}:I\to \mathbb{R}$$ be the elementwise constant function with $$\bar{q}|_{T_i}:=\sup_{x\in T_i}q(x)$$ for some non-negative function $$q:\bar{I}\to\mathbb{R}$$. The following bound holds: \begin{equation}\label{uhatmu_estimate2} \begin{aligned} | \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|_q^2 \le &\, \sum_{i\in \mathcal{N}_-} \left( C\left(\varepsilon h_i^{-1}+\bar{q}^2h_i\right)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) +\frac{|b|}{2}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) \right) \end{aligned} \end{equation} (3.4) for some constant $$C>0$$, independent of $$\varepsilon$$, $$b$$, $$c$$, $$q$$, $$h$$, $$\hat{u}_h$$ and of $$u_h$$, but dependent on $$p$$. Proof. This follows by Lemma 3.3 and properties $$\hat{u}_h(\partial_+I)=u_h(\partial_+I)$$ and $$\hat{u}_h(\partial_-I)=0$$. □ Remark 3.5. We note that the second term on the right-hand side of (3.4) can be removed from both the left- and right-hand sides of (3.4), without affecting the validity of this bound. 4. A posteriori error bound We present an a posteriori error bound for $$\varepsilon>0$$ whose proof remains valid at the singular limit $$\varepsilon=0$$ also; the singular limit will be discussed in Section 4.2. In the interest of clarity of presentation, we shall first consider the simpler case of nonvanishing reaction coefficient $$c$$ in (2.1) (Section 4.1), before moving to the general case of small or vanishing reaction in Section 4.3. A comparison with known a posteriori bounds from the literature is given in Section 6.4. To this end, we begin by considering an auxiliary boundary value problem. For $$\varepsilon>0$$, we set $$\tilde{u}: = u + g$$, where $$g:I\to \mathbb{R}$$ is a linear function on $$I$$ with $$g(\partial_-I) = 0$$ and $$g(\partial_+I) = u_h(\partial_+I)$$. Then, $$\tilde{u}\in H^1(I)$$ satisfies the boundary value problem \begin{equation}\label{pde_mod} -\varepsilon \tilde{u} ''+b \tilde{u}'+c\tilde{u} =f+bg'+cg\quad\text{in } I,\qquad \tilde{u}(\partial_-I) = 0, \qquad \tilde{u}(\partial_+I) = u_h(\partial_+I). \end{equation} (4.1) For the singular limit case $$\varepsilon=0$$, we consider only the inflow boundary condition $$\tilde{u}(\partial_-I)=0$$, and we set $$\tilde{u}:= u$$. Upon defining $$\rho :=\tilde{u}-\hat{u}_h$$, we observe that $$\rho\in H^1_0(I)$$ when $$\varepsilon>0$$, and $$\rho\in H^1(I)$$ when $$\varepsilon=0$$. 4.1 The case $$\varepsilon>0$$ and $$c>0$$ Let $$\pi:L^2(I)\to \mathbb{V}_p$$ denote the orthogonal $$L^2$$-projection onto the finite element space. Then, we have the following error equation for the reconstruction. Lemma 4.1 (Error equation for the reconstruction). Set $$z:=bg'+cg$$ for brevity. For $$\varepsilon>0$$, we have \begin{equation}\label{err_eq} \begin{aligned} \varepsilon (\rho',v') + B_{\rm c}(\rho,v) = (\,f-\pi f+z,v) - \varepsilon (\hat{u}'_h,v') +\varepsilon B_{\rm d}(u_h,\pi v) +(\pi(cu_h)-c\hat{u}_h,v) \end{aligned} \end{equation} (4.2) for all $$v\in H^1_0(I)$$. For $$\varepsilon =0$$, (4.2) holds for all $$v\in H^1(I)$$. Proof. From the continuity of the reconstruction $$\hat{u}_h$$ on $$I$$, we have \begin{equation}\label{err_one} \begin{aligned} \varepsilon (\rho',v') + (b\rho'+c\rho,v) &= (\,f+z,v) - \varepsilon (\hat{u}'_h,v')- (b\hat{u}'_h+c\hat{u}_h,v) \\ &= (\,f+z,v) - \varepsilon (\hat{u}'_h,v') - (b\hat{u}'_h,\pi v) -(c\hat{u}_h,v) , \end{aligned} \end{equation} (4.3) using (4.1) and the fact that $$\hat{u}'_h\in \mathbb{V}_p$$, respectively. Using (3.2) and (2.4), the third term on the right-hand side of (4.3) can be written as \begin{equation}\label{err_three} (b\hat{u}'_h,\pi v) = B_{\rm c}(u_h,\pi v) -(cu_h,\pi v) =-\varepsilon B_{\rm d}(u_h,\pi v) + (\,f,\pi v) -(cu_h,\pi v). \end{equation} (4.4) Applying (4.4) on (4.3), along with standard properties of the orthogonal $$L^2$$-projection, results in the right-hand side being equal to \[ (\,f-\pi f+z,v) - \varepsilon \left(\hat{u}'_h,v'\right) +\varepsilon B_{\rm d}(u_h,\pi v) +(\pi(cu_h)-c\hat{u}_h,v), \] which already implies the result. □ For the proof of the a posteriori bound for $$\varepsilon>0$$, the following extension of the bilinear form $$B_{\rm d}$$ from $$\mathbb{V}_p\times \mathbb{V}_p$$ to $$(H^1(I)+\mathbb{V}_p)\times ( H^1(I)+\mathbb{V}_p)$$, given by \[ \begin{aligned} B_{\rm d}(w,v) := & \sum_{i=1}^N (w',v')_{T_i}-\sum_{i=0}^{N} \left(\{ \kern -1.6mm \{{(\pi w)'}\} \kern -1.6mm \}[ \kern -.7mm [{v}] \kern -.7mm ]+ \{ \kern -1.6mm \{{(\pi v)'}\} \kern -1.6mm \}[ \kern -.7mm [{w}] \kern -.7mm ]-\sigma[ \kern -.7mm [{w}] \kern -.7mm ][ \kern -.7mm [{v}] \kern -.7mm ]\right)(x_i),\\ \end{aligned} \] will be useful. The extended $$B_{\rm d}$$ is both coercive and continuous in $$(H^1(I)+\mathbb{V}_p)\times ( H^1(I)+\mathbb{V}_p)$$ with respect to the $$| \kern -.25mm \|{\cdot}| \kern -.25mm \|$$-norm for sufficiently large penalty constant $$C_{\sigma}$$ (see, e.g., Georgoulis & Lasis, 2006 for a proof). Recalling the definition of $$\bar{q}:I\to \mathbb{R}$$ for a function $$q$$ from Lemma 3.4 (applied to $$c$$ below) we have the first of the main results of this work. Theorem 4.2. For $$\varepsilon>0$$ and $$c>0$$, the following a posteriori bound holds: \begin{align}\label{apost_final} | \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\sqrt{c}}^2& \le C\Bigg(\big\|{c_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\big\|_{}^2\nonumber\\ &\qquad\quad{} + \sum_{i=1}^{N-1} \!\! \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) +\sum_{i\in \mathcal{N}_-}\big((\varepsilon\sigma+ \bar{c} h)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2\big)(x_i)\nonumber\\ &\qquad\quad{} +((\varepsilon\sigma+ c_{\rm out} |b|)u_h^2)(\partial_+I)\Bigg) +\frac{|b|}{2}\sum_{i\in \mathcal{N}_-} [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i), \end{align} (4.5) with $$c_{\rm osc} := \min\{c^{-1/2}, \varepsilon^{-1/2}h\}$$, $$c_{\rm out}:= \max\{\|{cb^{-1}}\|_{},\|{bc^{-1}}\|_{}\}$$ and a constant $$C$$, independent of $$f$$, $$\varepsilon$$, $$b$$, $$c$$, $$h$$ and of $$u_h$$, but dependent on $$p$$ and on the size of the computational domain $$\beta-\alpha$$. Proof. Upon setting $$v=\rho$$ on (4.2), the terms on the left-hand side of (4.2) give \begin{equation} \varepsilon (\rho',v') + (b\rho'+c\rho,v) = \varepsilon \|{\rho'}\|_{}^2 + \|{\sqrt{c}\rho}\|_{}^2, \end{equation} (4.6) while the terms on the right-hand side of (4.2) can be bounded as follows. Since for $$\varepsilon>0$$, we have $$\hat{u}_h\in H^1(I)$$, with $$\hat{u}_h(\partial_-I) =0$$, $$\hat{u}_h(\partial_+I) =u_h(\partial_+I)$$ and $$\rho\in H^1_0(I)$$, we can deduce \[ (\hat{u}'_h,\rho') = B_{\rm d}(\hat{u}_h,\rho)+{\rm sign}(b) ((\pi \rho)'u_h)(\partial_+I), \] which implies, upon straightforward manipulation, \begin{align}\label{apost_one} B_{\rm d}(u_h,\pi\rho) -(\hat{u}'_h,\rho') & = B_{\rm d}(u_h,\pi\rho-\rho) + B_{\rm d}(u_h-\hat{u}_h,\rho)\nonumber\\ &\quad -{\rm sign}(b) ((\pi \rho)'u_h)(\partial_+I). \end{align} (4.7) For the first term on the right-hand side of (4.7), we set $$\tilde{v} =\pi\rho-\rho$$ (noting that $$\pi \tilde{v} =0$$ by construction) and we perform integration by parts on the elemental integrals: \[ \begin{aligned} B_{\rm d}(u_h,\tilde{v}) = & \sum_{i=1}^N (u'_h,\tilde{v}')_{T_i}-\sum_{i=0}^{N} \left(\{ \kern -1.6mm \{{u'_h}\} \kern -1.6mm \}[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]-\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ][ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]\right)(x_i)\\ = & \sum_{i=1}^{N-1} \left( [ \kern -.7mm [{u'_h\tilde{v}}] \kern -.7mm ]- \{ \kern -1.6mm \{{u'_h}\} \kern -1.6mm \}[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]+\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ][ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]\right)(x_i)+\sigma u_h\tilde{v}(x_N)+\sigma u_h\tilde{v}(x_0) , \end{aligned} \] since $$(u''_h,\tilde{v})_{T_i} =0$$ for all $$i=1,\dots,N$$, by the orthogonality of $$\pi$$. Further, noting the formula $$[ \kern -.7mm [{u'_h\tilde{v}}] \kern -.7mm ] = [ \kern -.7mm [{u'_h}] \kern -.7mm ]\{ \kern -1.6mm \{{\tilde{v}}\} \kern -1.6mm \}+\{ \kern -1.6mm \{{u'_h}\} \kern -1.6mm \}[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]$$ for all interior nodes, we arrive at \begin{equation}\label{diff_one} B_{\rm d}(u_h,\tilde{v}) = \sum_{i=1}^{N-1} \left([ \kern -.7mm [{u'_h}] \kern -.7mm ]\{ \kern -1.6mm \{{\tilde{v}}\} \kern -1.6mm \}+\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ][ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]\right)(x_i)+\sigma u_h\tilde{v}(x_N)+\sigma u_h\tilde{v}(x_0). \end{equation} (4.8) Using the approximation properties of $$\pi$$, we have \[ |[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ](x_i) |\le C \|{ \sqrt{h} \rho'}\|_{(x_{i-1},x_{i+1})},\quad |\{ \kern -1.6mm \{{\tilde{v}}\} \kern -1.6mm \}(x_i) |\le C \|{ \sqrt{h} \rho'}\|_{(x_{i-1},x_{i+1})}, \] and $$|\tilde{v}(x_0) |\le C \|{ \sqrt{h} \rho'}\|_{T_1}$$, $$|\tilde{v}(x_N) |\le C \|{ \sqrt{h} \rho'}\|_{T_N}$$, for some (generic) constant $$C>0$$, independent of $$\rho$$ and of $$h$$. Hence, the left-hand side of (4.8) can be estimated by \begin{equation}\label{diff_oneone} \begin{aligned} |B_{\rm d}(u_h,\tilde{v})| \le C \left(\sum_{i=1}^{N-1} \left(h|[ \kern -.7mm [{u'_h}] \kern -.7mm ]|^2\right)(x_i) +\sum_{i=0}^{N} \left(\sigma |[ \kern -.7mm [{u_h}] \kern -.7mm ]|^2\right)(x_i)\right)^{1/2}\|{\rho'}\|_{} \\ \end{aligned} \end{equation} (4.9) respectively, for some (generic) constant $$C>0$$, independent of $$\rho$$ and of $$h$$. The second term on the right-hand side of (4.7) can be estimated using the continuity of the extended bilinear form: \begin{equation}\label{diff_two} |B_{\rm d}(u_h-\hat{u}_h,\rho)|\le C | \kern -.25mm \|{u_h-\hat{u}_h}| \kern -.25mm \| | \kern -.25mm \|{\rho}| \kern -.25mm \|. \end{equation} (4.10) For the last term on the right-hand side of (4.7), using a standard inverse estimate along with the stability of the local $$L^2$$-projection in the $$H^1$$-seminorm, we deduce \begin{equation}\label{diff_threenew} | ((\pi \rho)'u_h)(\partial_+I)| \le C \|{\rho'}\|_{}(h^{-1/2}|u_h|)(\partial_+I). \end{equation} (4.11) Also, from the orthogonality of the $$L^2$$-projection, we have \begin{equation}\label{diff_three} \begin{aligned} &|(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| \\ &\quad{} = |((\,f-cu_h)-\pi (\,f-cu_h),\rho-\chi) -(c(\hat{u}_h-u_h),\rho)| \\ &\quad{} \le \|{(\,f-cu_h)-\pi (\,f-cu_h)}\|_{}\|{\rho-\chi}\|_{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\|{\sqrt{c}\rho}\|_{} \end{aligned} \end{equation} (4.12) for any $$\chi\in\mathbb{V}_p$$. The choice $$\chi = 0$$ results in \begin{align}\label{diff_threeone} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le \big(\|{c^{-1/2}((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}\nonumber\\ &\quad{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\big) \|{\sqrt{c}\rho}\|_{}, \end{align} (4.13) whereas the choice $$\chi = \pi \rho$$ results in \begin{align}\label{diff_threetwo} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le C \big(\|{h((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}\|{\rho'}\|_{}\nonumber\\ &\quad{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\|{\sqrt{c}\rho}\|_{}\big). \end{align} (4.14) Setting $$c_{\rm osc} = \min\{c^{-1/2}, \varepsilon^{-1/2}h\}$$, the last two bounds can be combined into \begin{align}\label{diff_threeall} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le C\big(\|{c_{osc}((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}\nonumber\\ &\quad{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\big){| \kern -.25mm \|{\rho}| \kern -.25mm \|}_{\sqrt{c}}. \end{align} (4.15) Finally, recalling that $$z=bg'+cg$$, we also have \begin{equation}\label{diff_fournew} |(z,\rho)| \le \|{c^{-1/2}z}\|_{}\|{\sqrt{c}\rho}\|_{} \le Cc_{\rm out} \sqrt{|b|}|u_h|(\partial_+I)\|{\sqrt{c}\rho}\|_{}. \end{equation} (4.16) Hence, using (4.9), (4.10), (4.11), (4.15) and (4.16) on (4.2) with $$v=\rho$$, we arrive at \begin{align}\label{apost_prefinal} {| \kern -.25mm \|{\rho}| \kern -.25mm \|}_{\sqrt{c}}^2 = \varepsilon \|{\rho'}\|_{}^2+\|{\sqrt{c}\rho}\|_{}^2 & \le C\Bigg(\|{c_{osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2\nonumber\\ &\qquad{} +\varepsilon| \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|^2+\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}^2\nonumber\\ &\qquad{} +\sum_{i=1}^{N-1} (\varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2)(x_i) +\sum_{i=0}^{N} (\varepsilon\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\nonumber\\ &\qquad{} + c_{\rm out} |b|u_h^2(\partial_+I)\Bigg). \end{align} (4.17) Also, we have \begin{equation}\label{tzi} \begin{aligned} | \kern -.25mm \|{g}| \kern -.25mm \|_{\sqrt{c}}^2 = &\ \varepsilon\|{g'}\|_{}^2+ \|{\sqrt{c}g}\|_{}^2+ \left(\left(\varepsilon\sigma +\frac{|b|}{2}\right) u_h^2\right)(\partial_+I)\\ \le &\ C((\varepsilon\sigma +c_{\rm out}|b|) u_h^2)(\partial_+I). \end{aligned} \end{equation} (4.18) Finally, observing the identity $$u-u_h = -g+\rho+\hat{u}_h-u_h,$$ the triangle inequality implies \[ | \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\sqrt{c}}\le | \kern -.25mm \|{g}| \kern -.25mm \|_{\sqrt{c}} +| \kern -.25mm \|{\rho}| \kern -.25mm \|_{\sqrt{c}}+| \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|_{\sqrt{c}}. \] The result readily follows by combining (4.17), (4.18) and the reconstruction stability estimates from Lemma 3.4. □ Theorem 4.2 above assumes that the reaction coefficient $$c$$ is positive. Moreover, the relative size between the convection and the reaction coefficients determines the magnitude of the constant $$c_{\rm out}$$ and, therefore, may adversely affect the constant of the a posteriori bound. It is possible to remove these restrictive effects, by considering a variant of Theorem 4.2; this will be the content of Section 4.3. 4.2 The singular limit case $$\varepsilon=0$$ Crucially, the bound (4.5) remains valid at the singular limit $$\varepsilon =0$$. An inspection of the proof of Theorem 4.2, noting carefully that in this case we have $$g=z=0$$, gives \[ \begin{aligned} & \|{\sqrt{c}(u-u_h)}\|_{}^2 + \frac{|b|}{2} \sum_{i\in\mathcal{N}_-}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) + \frac{|b|}{2} (u-u_h)^2(\partial_+I)\\ &\quad{} \le C\left( \|{c^{-1/2} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 +\sum_{i\in\mathcal{N}_-} (\bar{c} h[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i) \right) + \frac{|b|}{2} \sum_{i\in\mathcal{N}_-}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i). \end{aligned} \] In this case, we also trivially have an a posteriori bound for the $$L^2$$-norm of the error \begin{equation}\label{singular_limit_apost} \|{\sqrt{c}(u-u_h)}\|_{}^2 \le \, C\left( \|{c^{-1/2} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 +\sum_{i\in\mathcal{N}_-} (\bar{c} h[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\right)\!; \end{equation} (4.19) cf. Makridakis & Nochetto (2006), for the respective result in the context of dG time stepping. 4.3 The case of small or vanishing reaction $$c$$ Now, we shall remove the assumption on the magnitude of the reaction coefficient by considering a variant of Theorem 4.2, which will be shown using a modified test function. The use of such testing in a priori analysis of dG methods is classical (see, e.g., Johnson & Pitkäranta, 1986). Theorem 4.3. Let $$\varepsilon>0$$, and assume that $$\varepsilon$$, $$b$$, $$c$$ are such that $$c+|b|-\varepsilon>0$$. Further, we set $$\gamma:I\to\mathbb{R}_+$$ with $$\gamma:=(c+|b|-\varepsilon)^{1/2}$$. Then, the following a posteriori bound holds: \begin{align}\label{apost_final_tw} | \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}^2 & \le C\Bigg(\|{\tilde{c}_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2\nonumber \\ &\qquad{} + \sum_{i=1}^{N-1} \!\! \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) +\sum_{i\in \mathcal{N}_-}\big((\varepsilon\sigma+ \bar{\gamma}^2 h)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2\big)(x_i) \nonumber\\ &\qquad{} +((\varepsilon\sigma+ c_{\rm out,2} |b|)u_h^2)(\partial_+I)\Bigg) +\frac{|b|}{2}\sum_{i\in \mathcal{N}_-} [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i), \end{align} (4.20) with $$\tilde{c}_{\rm osc} := \min\{\gamma^{-1}, \varepsilon^{-1/2}h\}$$, $$c\rm _{out,2}:= \max\{\|{\gamma^2 b^{-1}}\|_{},\|{b\gamma^{-2}}\|_{}\}$$ and a constant $$C$$, independent of $$f$$, $$\varepsilon$$, $$b$$, $$c$$, $$h$$ and of $$u_h$$, but dependent on $$p$$ and on $$\beta-\alpha$$. Moreover, in the singular limit $$\varepsilon =0$$, we have \begin{equation}\label{singular_limit_apost_two} \|{\gamma(u-u_h)}\|_{}^2 \le \, C\left( \|{\gamma^{-1} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 +\sum_{i\in\mathcal{N}_-} (\bar{\gamma}^2 h[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\right)\!, \end{equation} (4.21) which is valid also in the case of the pure advection problem, i.e., when $$c=0$$. Proof. We begin by defining the weight function $$\vartheta(x):= e^{-{\rm sign}(b)x}$$, $$x\in I$$ and we observe the trivial identity $$\vartheta'=-{\rm sign}(b)\vartheta$$. Further, we define the $$\vartheta$$-weighted $$L^2$$-norm by $$\|{(\cdot)}\|_{\vartheta}:=(\int_I\vartheta(x)(\cdot)^2(x)\,\mathrm{d} x)^{1/2}$$. Upon setting $$v=\rho \vartheta$$, on (4.2), and noting that $$v'= \rho'\vartheta + \rho\theta'=\rho'\vartheta-{\rm sign}(b)\rho\vartheta$$, the terms on the left-hand side of (4.2) give \begin{align} \varepsilon (\rho',v') + (b\rho'+c\rho,v) & = \varepsilon \|{\rho'}\|_{\vartheta}^2 + \|{\sqrt{c}\rho}\|_{\vartheta}^2 +(b-{\rm sign}(b) \varepsilon)(\rho',\rho\vartheta)\nonumber\\ & = \varepsilon \|{\rho'}\|_{\vartheta}^2 + \|{\gamma \rho}\|_{\vartheta}^2, \end{align} (4.22) upon integration by parts. The terms on the right-hand side of (4.2) can be bounded in the spirit of the proof of Theorem 4.2; we shall briefly repeat the estimates below focusing on the differences. Noting that $$v=\rho\vartheta\in H^1_0(I)$$, we have, as before, \begin{align}\label{apost_one_two} B_{\rm d}(u_h,\pi v) -(\hat{u}'_h,v') & = B_{\rm d}(u_h,\pi v-v) + B_{\rm d}(u_h-\hat{u}_h,v)\nonumber\\ &\quad{} -{\rm sign}(b) ((\pi v)'u_h)(\partial_+I). \end{align} (4.23) The first term on the right-hand side of (4.23), setting $$\tilde{v} =\pi v-v$$, performing integration by parts and using the approximation properties of $$\pi$$, as before, can be estimated by \begin{equation}\label{diff_oneone_two} |B_{\rm d}(u_h,\tilde{v})| \le C \left(\sum_{i=1}^{N-1} (h|[ \kern -.7mm [{u'_h}] \kern -.7mm ]|^2)(x_i) +\sum_{i=0}^{N} (\sigma |[ \kern -.7mm [{u_h}] \kern -.7mm ]|^2)(x_i)\right)^{1/2}\|{\rho'}\|_{\vartheta}, \end{equation} (4.24) observing the straightforward bound $$\|{v'}\|_{}\le C (\|{\rho'}\|_{\vartheta}+\|{\rho}\|_{\vartheta})\le C \|{\rho'}\|_{\vartheta}$$, with $$C$$ depending only on the interval $$I$$. The second term on the right-hand side of (4.23) can be estimated using the continuity of the extended bilinear form: \begin{equation}\label{diff_two_two} |B_{\rm d}(u_h-\hat{u}_h,v)|\le C | \kern -.25mm \|{u_h-\hat{u}_h}| \kern -.25mm \|\|{\rho'}\|_{\vartheta} , \end{equation} (4.25) since $$| \kern -.25mm \|{v}| \kern -.25mm \|_d=\|{v'}\|_{}\le C \|{\rho'}\|_{\vartheta}$$. For the last term on the right-hand side of (4.23), using a standard inverse estimate, we deduce \begin{equation}\label{diff_threenew_two} |((\pi v)'u_h)(\partial_+I)| \le C \|{v'}\|_{}(h^{-1/2}|u_h|)(\partial_+I)\le C \|{\rho'}\|_{\vartheta}(h^{-1/2}|u_h|)(\partial_+I). \end{equation} (4.26) Also, from the orthogonality of the $$L^2$$-projection, and working as before, we can have \begin{align}\label{diff_threeall_two} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le C\big(\|{\tilde{c}_{\rm osc}((\,f-cu_h)-\pi (\,f-cu_h))}\|_{\vartheta}\nonumber\\ &\quad{} +\|{\gamma(\hat{u}_h-u_h)}\|_{\vartheta}\big)( \varepsilon \|{\rho'}\|_{\vartheta}^2+ \|{\gamma \rho}\|_{\vartheta}^2)^{1/2}, \end{align} (4.27) with $$\tilde{c}_{osc} = \min\{\gamma^{-1/2}, \varepsilon^{-1/2}h\}$$. Finally, we also have \begin{equation}\label{diff_fournew_two} |(z,v)| \le \|{\gamma^{-1}z}\|_{\vartheta}\|{\gamma \rho}\|_{\vartheta} \le Cc_{\rm out,2}\sqrt{|b|}|u_h|(\partial_+I)\|{\gamma\rho}\|_{\vartheta}. \end{equation} (4.28) Using the above bounds, we can arrive at \begin{align}\label{apost_prefinal_two} \varepsilon \|{\rho'}\|_{\vartheta}^2+\|{\gamma\rho}\|_{\vartheta}^2 & \le C\Bigg(\|{\tilde{c}_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2\nonumber\\ &\quad{} +\varepsilon| \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|^2+\|{\gamma(\hat{u}_h-u_h)}\|_{}^2\nonumber\\ &\quad{} +\sum_{i=1}^{N-1} (\varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2)(x_i) +\sum_{i=0}^{N} (\varepsilon\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\nonumber\\ &\quad{} +c_{\rm out,2} |b|u_h^2(\partial_+I)\Bigg). \end{align} (4.29) Also, analogously to (4.18), we have $$| \kern -.25mm \|{g}| \kern -.25mm \|_{\gamma}^2 \le C((\varepsilon\sigma +c_{\rm out,2}|b|) u_h^2)(\partial_+I)$$. Finally, as before, the triangle inequality and the reconstruction stability estimates from Lemma 3.4 provide the first bound. The bound in the singular limit $$\varepsilon =0$$ follows by an inspection of the proof above with $$\varepsilon=g=z=0$$. □ Remark 4.4. The above result improves upon Theorem 4.2, in terms of the dependence with respect to the relative size between the reaction and the convection coefficients. Indeed, an $$L^2$$-component remains positive when $$c=0$$ and the constant $$c\rm _{out,2}$$ does not degenerate as $$c\to 0$$. 4.4 On extension to two space dimensions Here, we comment on the current obstacles on extending the above analysis to two (or more) space dimensions. First and foremost, the results for the respective hyperbolic problem presented in Georgoulis et al. (2014) are quite preliminary at this point and are valid only for special structured meshes in sufficient generality for the above analysis to carry through. A second challenge, of a rather technical nature, is the extension of the construction (4.1) in two dimensions. 5. On lower bounds The above a posteriori bound (4.5) can be complemented by the following lower bound, Theorem 5.1. This is derived by utilizing standard techniques and, as expected, does not provide enough information for very small values of $$\varepsilon$$ and moderate values of $$h.$$ A more detailed analysis using the structure of the possible boundary layer yields the second result in this section, Theorem 5.2, which establishes the quality of the estimator for small values of $$\varepsilon$$. Theorem 5.1. For each element $$T_i\in\mathcal{T}$$, $$i=1,\dots,N$$, we have \begin{align}\label{second_lower} \varepsilon^2 h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & \le C\left(h^2\tilde{{\rm Osc}}_i^2 + \varepsilon\|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + |b|^2 \|{u-u_h}\|_{T_i\cup T_{i+1}}^2 \right.\nonumber \\ &\quad{} \left.+|b|^2 \left( h_i [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_{i-1})+h_{i+1} [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_{i})\right) + c_i \|{h\sqrt{c}(u-u_h)}\|_{T_i\cup T_{i+1}}^2\right)\!, \end{align} (5.1) with $${\rm Osc}_i: = \|{f-cu_h-\pi (\,f-cu_h)}\|_{T_i}$$, $$i=1,\dots,N$$ and $$\tilde{{\rm Osc}}_i :=({\rm Osc}_i^2+{\rm Osc}_{i+1}^2)^{1/2}$$, and $$C>0$$ a constant dependent only on $$p$$ and on the local quasi-uniformity of the mesh around $$x_i$$. Proof. We start by defining an elemental bubble function $$\beta_i$$ to be a quadratic polynomial on the element $$T_i$$, such that $$\beta_i(x_{i-1}) = 0 = \beta_i(x_i)$$ and $$\beta_i((x_{i-1}+x_i)/2) = 1$$, resulting in $$0\le \beta_i\le 1$$ and $$\|{\beta_i}\|_{T_i}^2 = C h_i$$, for a constant $$C$$ independent of $$h_i$$. For $$\psi : = h_i^2\phi\beta_i$$, with $$\phi:= \pi f +\varepsilon u_h'' -b\hat{u}_h'-\pi(c u_h)$$, we have, respectively, \begin{align} & \int_{T_i} \varepsilon (u-u_h)'\psi' -b(u-\hat{u}_h)\psi'+c(u-u_h)\psi \,\mathrm{d} x\nonumber\\ &\quad{} = \int_{T_i}\phi\psi \,\mathrm{d} x - \int_{T_i}(\,f-cu_h-\pi(\,f-cu_h))\psi \,\mathrm{d} x, \end{align} (5.2) which gives \begin{align} h_i^2\|{\sqrt{\beta_i} \phi}\|_{T_i}^2 & \le h_i^2 {\rm Osc}_i\|{\phi}\|_{T_i} + C \varepsilon h_i \|{(u-u_h)'}\|_{T_i}\|{\phi}\|_{T_i}\nonumber\\ &\quad{} + C |b| h_i\|{u-\hat{u}_h}\|_{T_i}\|{\phi}\|_{T_i}\nonumber\\ &\quad{} + \sqrt{c_i} h_i^2 \|{\sqrt{c}(u-u_h)}\|_{T_i}\|{\phi}\|_{T_i}, \end{align} (5.3) and, finally, using the (standard) bound $$\|{\phi}\|_{T_i} \le C \|{\sqrt{\beta_i} \phi}\|_{T_i}$$, for some $$C$$ independent of $$\phi$$ and of $$T_i$$, we get \begin{equation}\label{first_lower} \begin{aligned} h_i^2\|{ \phi}\|_{T_i}^2 \le &\ C\left( h_i^2 {\rm Osc}_i^2+ \varepsilon \|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i}^2 \right.\\ &\left.+ |b|^2\|{u-\hat{u}_h}\|_{T_i}^2 + c_i h_i^2 \|{\sqrt{c}(u-u_h)}\|_{T_i}^2\right)\!. \end{aligned} \end{equation} (5.4) Having estimated $$\|{ \phi}\|_{T_i}$$ by a norm of the error, we shall now estimate the term $$\varepsilon[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2$$ from above by the same norm of the error. To this end, we define the bubble function $$\tilde{\beta}_i$$ to be a linear polynomial on each $$T_i$$ and $$T_{i+1}$$, such that $$\tilde{\beta}_i(x_{i-1}) = 0 = \tilde{\beta}_i(x_{i+1})$$ and $$\tilde{\beta}_i(x_i) = 1$$ on $$T_i$$, resulting in $$0\le \tilde{\beta}_i\le 1$$ and $$\|{\tilde{\beta}_i}\|_{T_j}^2 = h_j/2$$, $$j=i,i+1$$. Now, for $$\tilde{\psi} : = \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\tilde{\beta}_i$$, we have, respectively, \begin{align} & \int_{T_i\cup T_{i+1}} \!\!\varepsilon (u-u_h)'\tilde{\psi}' -b(u-\hat{u}_h)\tilde{\psi}'+c(u-u_h)\tilde{\psi} \,\mathrm{d} x\nonumber \\ &\quad{} = \int_{T_i\cup T_{i+1}}\!\!\! \phi\tilde{\psi} \,\mathrm{d} x - \varepsilon([ \kern -.7mm [{u'_h}] \kern -.7mm ]\tilde{\psi})(x_i), \end{align} (5.5) noting that $$\int_{T_i\cup T_{i+1}}\!(\,f-cu_h-\pi(\,f-cu_h))\tilde{\psi} \,\mathrm{d} x =0$$, because $$\tilde{\psi} \in \mathbb{V}_p$$. This gives \begin{align} \varepsilon^2 h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & \le \|{\phi}\|_{T_i\cup T_{i+1}}\varepsilon h^{3/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\nonumber\\ &\quad{} + C\|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i\cup T_{i+1}}\varepsilon^{3/2} h^{1/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\nonumber\\ &\quad{} + C |b| \|{u-\hat{u}_h}\|_{T_i\cup T_{i+1}}\varepsilon h^{1/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\nonumber\\ &\quad{} + \sqrt{c_i} \|{\sqrt{c}(u-u_h)}\|_{T_i\cup T_{i+1}}\varepsilon h^{3/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i), \end{align} (5.6) and, thus, \begin{align}\label{second_lower2} \varepsilon^2 h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & \le C\left( \|{h\phi}\|_{T_i\cup T_{i+1}}^2 + \varepsilon\|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 \right.\nonumber\\ &\quad{} \left.+ |b|^2 \|{u-\hat{u}_h}\|_{T_i\cup T_{i+1}}^2 + c_i \|{h\sqrt{c}(u-u_h)}\|_{T_i\cup T_{i+1}}^2\right)\!. \end{align} (5.7) Using (5.4) to estimate the first term on the right-hand side of (5.7), and using Lemma 3.3, we finally arrive at the result. □ We remark that the remaining terms in the a posteriori bound above either involve jump-type terms of the numerical solution and are present also in the dG-energy norm or represent data oscillation. Indeed, it is possible to add the remaining jump-type terms from Theorems 4.2 and 4.3 to both sides of (5.1) to arrive at a complete lower bound. We note that the constant in the above lower bound is proportional to $$\varepsilon$$. This is to be expected as the continuity constant of the bilinear form $$\epsilon B_{\rm d} +B_{\rm c}$$ with respect to the dG-energy norm $$| \kern -.25mm \|{\cdot}| \kern -.25mm \|_{\gamma}$$ depends unfavourably on the Péclet number. On the other hand, as we shall see in the numerical experiments below, the term $$\varepsilon h [ \kern -.7mm [{u_h'}] \kern -.7mm ]^2$$ does not appear to be dominant in the pre-asymptotic regime (i.e., when $$h> \varepsilon$$) at least for problems involving boundary and interior layers. This observation can be further substantiated as follows. Let $$u_0\in H^1(I)$$ be the exact solution to the reduced problem $$bu'_0+cu_0 = f$$, $$u_0(\partial_- I)=0$$, and define $$s_0\in C^1(\bar{I})$$ to be a $$C^1$$-piecewise polynomial interpolant of $$u_0$$ (e.g., Hermite interpolant) with nodes $$x_i$$, $$i=0,1\dots,N$$. Also, let $$\pi_0:L^2(I)\to \mathbb{V}_0$$ denote the orthogonal $$L^2$$-projection operator onto constants on each $$T_i$$. We then have \begin{align}\label{lower_bound_two} \varepsilon h_i[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & =\varepsilon h_i[ \kern -.7mm [{(s_0-u_h)'}] \kern -.7mm ]^2(x_i)\nonumber\\ &\le C\varepsilon \sum_{j=i,i+1}\|{(s_0-u_h)'}\|_{T_j}^2\nonumber\\ &=C\varepsilon \sum_{j=i,i+1}\|{(s_0-u_h-\pi_0(u-u_h))'}\|_{T_j}^2\nonumber\\ &\le C\varepsilon \sum_{j=i,i+1}h_j^{-2}\|{s_0-u_h-\pi_0(u-u_h)}\|_{T_j}^2\nonumber\\ &\le C\varepsilon \sum_{j=i,i+1}h_j^{-2}\Big( \|{u-u_h-\pi_0(u-u_h)}\|_{T_j}^2+ \|{u-u_0}\|_{T_j}^2+ \|{u_0-s_0}\|_{T_j}^2\Big)\nonumber\\ &\le C \sum_{j=i,i+1}\Big( \varepsilon \|{(u-u_h)'}\|_{T_j}^2+ \frac{\varepsilon }{h_j}\|{u-u_0}\|_{L^{\infty}(T_j)}^2+ \frac{\varepsilon }{h_j^{2}}\|{u_0-s_0}\|_{T_j}^2\Big), \end{align} (5.8) using inverse estimates and the local quasi-uniformity of the mesh, respectively. Now, classical asymptotic expansions (see, e.g., Roos et al., 2008, Theorem I.1.4) ensure that there exists a constant $$C>0$$, independent of $$\varepsilon$$, such that \begin{equation}\label{asymptotic} \|{u-u_0}\|_{L^{\infty}(I_{\zeta})}\le C\varepsilon, \end{equation} (5.9) with $$I_{\zeta}:=\{x\in I: {\rm dist} (x, \partial_+I)>\zeta\}$$, for fixed $$\zeta>0$$, sufficiently large, so that it does not include the boundary layer region of size $$\mathcal{O}(\varepsilon)$$ in the vicinity of $$\partial_+I$$. Also, since $$\|u\|_{L^{\infty}(I)}\le C$$ uniformly with respect to $$\varepsilon$$ (see, e.g., Roos et al., 2008, Theorem I.1.4), we have \[ \|{u-u_0}\|_{T_j\cap (I\backslash I_\zeta)}^2\le C|T_j\cap (I\backslash I_\zeta)|\le Ch_j \] for all $$T_j$$ with $$T_j\cap (I\backslash I_\zeta)\ne \emptyset$$, with $$C>0$$ a constant independent of $$\varepsilon$$ and of $$h_j$$. (Note that taking $$\varepsilon$$ small enough, we can have only one element $$T_j$$ such that $$T_j\cap (I\backslash I_\zeta)\ne \emptyset$$ in the mesh.) Moreover, since $$u_0\in H^1(I)$$, with $$u_0(x_i)=s_0(x_i)$$ from the interpolation property of $$s_0$$ at the nodes, the Poincaré–Friedrichs inequality implies \begin{equation}\label{pf_reduced} \|{u_0-s_0}\|_{T_i}\le \frac{h_i}{\pi} \|{(u_0-s_0)'}\|_{T_i} \end{equation} (5.10) for $$i=1,\dots,N$$. Furthermore, selecting $$s_0$$ to be exactly the Hermite interpolant of $$u_0$$ at the nodes $$x_i$$ and, assuming also $$f,c\in H^1(T_i)$$, we have $$u_0\in H^2(T_i)$$ and thus, \begin{equation}\label{pf_reduced_two} \|{u_0-s_0}\|_{T_i}\le \frac{h_i^2}{\pi^2} \|{(u_0-s_0)''}\|_{T_i}\le C h_i^2, \end{equation} (5.11) for some constant $$C>0$$ (which is, by construction, independent of $$\varepsilon$$). Combining the above estimates, we have the following result. Theorem 5.2. Let $$u_0\in H^{1+r}(I)$$, for $$r= 0,1$$. Then, for each element $$T_i\in\mathcal{T}$$, $$i=1,\dots,N$$ with $$T_i\in I_\zeta$$, we have \begin{equation}\label{lower_two_thm} \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon^3 h_i^{-1} + \varepsilon h_i^{2r} \big), \end{equation} (5.12) with the constant $$C>0$$, independent of $$\varepsilon$$, $$h$$ and of $$u$$, but possibly dependent on $$b$$, $$c$$ and on $$f$$. Moreover, we have \begin{equation}\label{lower_two_thm_two} \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon h_i^{-1} + \varepsilon h_i^{2r} \big) \end{equation} (5.13) for $$(T_i\cup T_{i+1}) \cap (I\backslash I_\zeta)\ne \emptyset$$, also. □ Now, assume for the moment that $$h_i>\varepsilon$$ for all $$i=1,\dots,N$$, and let $$k$$ be the index of the element $$T_k$$ having the outflow boundary $$\partial_+I$$ as one of its end points (i.e., $$k\in\{1,N\}$$) and that $$\varepsilon$$ is small enough, so that choosing $$\zeta =h_k$$ ensures the validity of (5.9) for this $$I_\zeta$$. In this setting, Theorem 5.2 implies \[ \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon h_i^r \big) \] for $$T_i\cup T_{i+1} \subset I_\zeta$$, and \begin{equation}\label{boundary_layer_element} \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon h_i^{-1} \big) \end{equation} (5.14) for $$(T_i\cup T_{i+1}) \cap (I\backslash I_\zeta)\ne \emptyset$$, with $$k\in \{i,i+1\}$$. Therefore, fixing $$h$$ and letting $$\varepsilon\to 0$$ highlights the robustness of the proposed a posteriori estimator in the pre-asymptotic regime $$h_k>\varepsilon$$ for all elements in $$I_\zeta$$. Moreover, for $$\varepsilon$$ small enough, we can safely assume that $$T_k=I\backslash I_\zeta$$. For this element, the lower bound (5.14) shows that the estimator remains bounded as $$h_k\to \varepsilon$$. On the other hand, the good quality of the proposed estimator follows from Theorem 5.1 in the asymptotic regime $$h_k\le \varepsilon$$. In this latter case, dividing both sides of (5.1) by $$\varepsilon$$ and noticing that the last three terms on the right-hand side of (5.1) are of higher order with respect to the local mesh size, the boundedness of the local mesh-Péclet number $$|b|h_j/\varepsilon$$, $$j=i,i+1$$ ensures the good behaviour of the constant in (5.1). 6. Numerical experiments We present a number of examples to verify the quality of the estimators proposed in Theorems 4.2 and 4.3, especially their stability as the Péclet number $$|b|/\varepsilon\to \infty$$ and their ability to drive automatic adaptive mesh-refinement strategies. In Section 6.3, we discuss the a posteriori error indicator presented above, in conjunction with the augmented-norm-type estimators and in particular with the ones from Schötzau & Zhu (2009) and Ern et al. (2010). 6.1 Example 1 We begin by studying how the new a posteriori error estimate behaves for a problem with a smooth solution for various values of $$\varepsilon$$ under uniform mesh refinement. In this case, we choose $$(\alpha,\beta) = (0,1)$$, $$b=1$$, $$c=1$$ and pick $$f$$ such that $$u (x) = \sin(8\pi x)$$. Figure 1(a) shows the convergence of the error measured in the dG norm as the mesh is refined, while Fig. 1(b) shows the effectivities $$\eta/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ as the mesh is refined, for $$\varepsilon = 1,1e-1,\ldots,1e-4$$ and for $$\varepsilon = 0$$, with $$\eta$$ denoting the right-hand side of (4.20) without the constant, namely, \[ \begin{aligned} \eta:=&\|{\tilde{c}_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 + \sum_{i=1}^{N-1} \!\! \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \\ &+\sum_{i\in \mathcal{N}_-}\left(\left(\varepsilon\sigma+ \bar{\gamma}^2 h+\frac{|b|}{2}\right)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2\right)(x_i) +((\varepsilon\sigma+ c_{\rm out,2} |b|)u_h^2)(\partial_+I). \end{aligned} \] Fig. 1. View largeDownload slide Example 1. Smooth problem. (a) Error convergence under uniform refinement. (b) Effectivity of estimator under uniform refinement. Fig. 1. View largeDownload slide Example 1. Smooth problem. (a) Error convergence under uniform refinement. (b) Effectivity of estimator under uniform refinement. As expected from a priori error analysis (see, e.g., Houston et al., 2002) and noting that, for a uniform mesh, $$h \sim N^{-1}$$, Fig. 1(a) shows that for $$\varepsilon > 0$$, convergence is $$\mathcal{O}(h)$$ in the asymptotic regime, whereas for $$\varepsilon = 0$$, the order is $$\mathcal{O}(h^{3/2})$$. We note that, for small $$\varepsilon$$, the jump terms in the dG-energy norm dominate on coarse meshes, hence initially $$\mathcal{O}(h^{3/2})$$ convergence is witnessed. For all $$\varepsilon$$, the effectivities remain bounded between $$0.9$$ and $$3.5$$, showing robustness of the error estimator as $$\varepsilon \rightarrow 0$$. In fact, as the mesh is refined, effectivities tend to a little over $$3$$ for $$\varepsilon >0$$, whereas, effectivities tend to $$1$$ for $$\varepsilon = 0$$. 6.2 Example 2 In this second example, we again select $$b=1$$ and $$c=1$$, but let $$f=1$$ on the domain $$(\alpha,\beta) = (0,1)$$; the exact solution exhibits a boundary layer near $$x=1$$, which narrows as $$\varepsilon$$ is reduced. Once again, we test the robustness of the estimator $$\eta$$, i.e., the right-hand side of (4.20), but this time we do so within an adaptive strategy. Elements are marked for refinement using a bulk criterion, where those elements contributing to 50% of the total error are refined; no elements are chosen for coarsening. In all cases, the computation starts from a uniform mesh comprising $$8$$ elements. We test with $$\varepsilon = 1, 1e-1,\ldots,1e-7$$ and $$\varepsilon = 0$$, noting there is no boundary layer in the latter case. Figure 2(a) shows the effectivities $$\eta/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ as the mesh is adaptively refined. In all cases, effectivities remain bounded between $$1$$ and $$3.5$$ and converge to a value just over $$3$$ for $$\varepsilon > 0$$ and to $$1$$ for $$\varepsilon = 0$$, highlighting the robustness of the proposed estimator. Fig. 2. View largeDownload slide Example 2. Boundary layer problem under adaptive refinement. (a) Effectivity of the error estimator. (b) Error convergence under uniform and adaptive refinement. Fig. 2. View largeDownload slide Example 2. Boundary layer problem under adaptive refinement. (a) Effectivity of the error estimator. (b) Error convergence under uniform and adaptive refinement. Figure 2(b) compares the error convergence for the adaptive strategy against uniform refinement for $$\varepsilon = 1e-2$$ and $$\varepsilon = 1e-5$$. For $$\varepsilon = 1e-2$$, on the latter meshes, the adaptive strategy gives well over an order of magnitude improvement in error over uniform refinement for a comparable number of degrees of freedom. For $$\varepsilon = 1e-5$$ and for around 5,000 degrees of freedom, the uniform strategy has failed to resolve the boundary layer successfully, while the adaptive strategy is approximately $$3$$ orders of magnitude superior for comparable degrees of freedom. 6.3 Example 3 We now test the bound from Theorem 4.3 by selecting $$c=0$$. As before, we let $$(\alpha,\beta) = (0,1)$$, $$b=1$$ but choose $$f = {\rm exp}(-10,000(x-0.5)^2)$$, which introduces an interior layer into the solution, narrowing as $$\varepsilon$$ is reduced. To prevent a boundary layer forming, we choose a homogeneous Neumann boundary condition on the right-hand side of the domain $$u^\prime(1) = 0$$; the estimator from Theorem 4.3 is modified in a trivial fashion to allow for this. Once again, we start with a uniform mesh of $$8$$ elements and carry out refinement using the same bulk criterion marking strategy as in Example 2, in this case for $$\varepsilon = 1e-1,\dots,1e-7$$. Figure 3(a) shows the effectivities as the mesh is refined. We see that the effectivities remain bounded between $$1$$ and $$20$$. This is a larger interval than in both Examples 1 and 2, although the largest effectivity occurs for $$1e-1$$. Fig. 3. View largeDownload slide Example 3. Interior layer problem under adaptive refinement. (a) Effectivities for varying $$\varepsilon$$. (b) Error convergence under uniform and adaptive refinement. Fig. 3. View largeDownload slide Example 3. Interior layer problem under adaptive refinement. (a) Effectivities for varying $$\varepsilon$$. (b) Error convergence under uniform and adaptive refinement. A comparison between uniform and adaptive refinement is shown in Fig. 3(b) for $$\varepsilon = 1e-2$$ and $$\varepsilon = 1e-4$$. In both cases, on later meshes, the adaptive refinement algorithm shows around 1 order of magnitude improvement in error for comparable degrees of freedom. 6.4 Comparison with other a posteriori bounds A comparison with other a posteriori estimates is in place. The pioneering works Sangalli (2001), Verfürth (2005), Sangalli (2008), Schötzau & Zhu (2009) introduce augmented norms for which a posteriori estimators are shown to admit so-called robust upper and lower bounds. In the context of dG methods, the augmented norm \[ | \kern -.25mm \|{w}| \kern -.25mm \|^2_{\rm aug} = \varepsilon| \kern -.25mm \|{w}| \kern -.25mm \|^2+\|cw\|^2 +\sum_{i=0}^N\left(ch+\frac{h}{\varepsilon}\right)[ \kern -.7mm [{w}] \kern -.7mm ]^2(x_i)+|bw|_{\ast}^2, \] where $$|q|_{\ast} := \sup_{v \in H^1_0(\alpha,\beta)\backslash\{0\}}\big(\int_{\alpha}^{\beta} qv' \,\mathrm{d} x)/(\varepsilon| \kern -.25mm \|{v}| \kern -.25mm \|^2+\|cv\|^2\big)^{1/2}$$ for $$q \in L^2(\alpha,\beta)$$ was used in Schötzau & Zhu (2009) to show that the estimator \begin{align}\label{sz} \eta^2_{\rm SZ} & = \sum_{i=1}^N\|c_{\rm osc}(\pi f+\varepsilon u_h^{\prime\prime}-bu^\prime_h-cu_h)\|^2_{T_i}+\sum_{i=1}^{N-1}\varepsilon^{-1/2}c_{\rm osc}[ \kern -.7mm [{u^\prime_h}] \kern -.7mm ]^2(x_i)\nonumber\\ &\quad{} +\sum_{i=0}^N\left(\varepsilon \sigma+ch+\frac{h}{\varepsilon}\right)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) \end{align} (6.1) is robust for $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\rm aug}^2$$, up to data oscillation; here we have taken $$c\ge 0$$ constant, for simplicity of the exposition. The inclusion of a dual-type $$|\cdot|_*$$-seminorm is an essential feature for the robustness of the estimator (6.1); corresponding observations are also true in the context of stabilized conforming methods (Sangalli, 2001, 2008; Verfürth, 2005). Indeed, going back to the setting of Example 1, the respective effectivities $${\eta}_{\rm SZ}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ with $${\eta}_{\rm SZ}$$ as in (6.1) and $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ the dG-energy norm defined in (2.5) (which, crucially, does not include the $$|\cdot|_*$$-seminorm) grow as $$\varepsilon$$ becomes smaller in the pre-asymptotic regime. In Fig. 4(a), we give these effectivities against the number of degrees of freedom $$N$$. Figures 4(b,c) show $$\eta_{\rm R}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ and $$\eta_{\rm J}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ against $$N$$, respectively, with \[ \eta^2_{\rm R}=\sum_{i=1}^Nc_{\rm osc}\|\pi f+\varepsilon u_h^{\prime\prime}-bu^\prime_h-cu_h\|^2_{T_i}\qquad\text{ and }\qquad \eta^2_{\rm J}=\sum_{i=0}^{N}\frac{h}{\varepsilon}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i). \] Fig. 4. View largeDownload slide Example 1. Effectivity of the a posteriori estimator (6.1) under uniform refinement without the inclusion of the $$|\cdot|_*$$ in the error norm. (a) Effectivity of the $$\eta_{SZ}$$ estimator. (b) Contribution of the residual term $$\eta_{\rm R}$$. (c) Contribution of the jump term $$\eta_{\rm J}$$. Fig. 4. View largeDownload slide Example 1. Effectivity of the a posteriori estimator (6.1) under uniform refinement without the inclusion of the $$|\cdot|_*$$ in the error norm. (a) Effectivity of the $$\eta_{SZ}$$ estimator. (b) Contribution of the residual term $$\eta_{\rm R}$$. (c) Contribution of the jump term $$\eta_{\rm J}$$. Crucially, the estimator from Theorems 4.2 and 4.3 does not contain the terms $$\eta_{\rm R}$$ and $$\eta_{\rm J}$$ and, hence, robustness with respect to the standard dG-energy norm $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ is possible. To assess the potential impact of the above discussion in the context of adaptivity, we compare using $$\eta_{\rm SZ}$$ or $$\eta$$ to drive the adaptive algorithm for the problem of Example 2. A comparison of using $$\eta_{\rm SZ}$$ or $$\eta$$ to drive the adaptive algorithm is carried out for various values of $$\varepsilon$$ with the error measured in the dG-energy norm $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$; in Fig. 5, we show the results for $$\varepsilon=1e-5,1e-6,1e-7$$. We notice that $$\eta$$ initially leads to meshes that give a reduced error for a comparable number of degrees of freedom; indeed, the improvement is more marked for smaller $$\varepsilon$$. Fig. 5. View largeDownload slide Example 2. Performance of adaptive algorithm driven by $$\eta$$ and by $$\eta_{\rm SZ}$$. (a) $$\varepsilon = 1e-5$$. (b) $$\varepsilon = 1e-6$$. (c) $$\varepsilon = 1e-7$$. Fig. 5. View largeDownload slide Example 2. Performance of adaptive algorithm driven by $$\eta$$ and by $$\eta_{\rm SZ}$$. (a) $$\varepsilon = 1e-5$$. (b) $$\varepsilon = 1e-6$$. (c) $$\varepsilon = 1e-7$$. We also consider the estimator based on reconstructions of the diffusive and convective fluxes onto the Raviart–Thomas finite element space and on reconstructions of the potential into $$H^1_0(\it\Omega)$$, presented in Ern et al. (2010). Let the restriction of this estimator to the one-dimensional setting be denoted by $$\eta_{\rm ESV}$$; we omit full details of $$\eta_{\rm ESV}$$ for conciseness and we refer to Ern et al. (2010, Theorem 3.5). We investigate numerically whether $$\eta_{\rm ESV}$$ is also a robust estimator for the dG-energy norm in the present one-dimensional setting for the problem of Example 1. Figure 6 shows the effectivities $${\eta}_{\rm ESV}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ for $$\epsilon = 1\times 10^n$$, $$n=-4,\ldots,0$$, as the mesh is uniformly refined. We notice that, as $$\epsilon$$ is reduced, the maximum effectivity increases as in the case of $$\eta_{SZ}$$. Fig. 6. View largeDownload slide Effectivity of estimator under uniform refinement. (a) Effectivity of $$\eta_{\rm ESV}$$ under uniform refinement. (b) Effectivity of the element residual term $$\eta_{R,T}$$ (cf. Ern et al., 2010, Equation (35) for precise definition) of $$\eta_{\rm ESV}$$ under uniform refinement. Fig. 6. View largeDownload slide Effectivity of estimator under uniform refinement. (a) Effectivity of $$\eta_{\rm ESV}$$ under uniform refinement. (b) Effectivity of the element residual term $$\eta_{R,T}$$ (cf. Ern et al., 2010, Equation (35) for precise definition) of $$\eta_{\rm ESV}$$ under uniform refinement. 6.5 Example 4 To highlight the potential generality of the proposed estimator for higher-dimensional problems, we consider an error estimate analogous to that of Theorem 4.2 in the two-dimensional setting. As mentioned in the introduction, such a result can be obtained, provided we have at our disposal a technique that leads to error control of the limiting hyperbolic problem. Such results are available in certain cases (Georgoulis et al., 2014). To this end, we consider the problem \begin{eqnarray} -\varepsilon{\it\Delta} u+{\bf b}\cdot \nabla u + u &=& 1, \quad (x,y) \in {\it\Omega} = (0,1)^2, \nonumber\\ u &=& 0, \quad (x,y) \in \partial {\it\Omega}, \nonumber \end{eqnarray} where $${\bf b} = (\sqrt{2},1)^{\rm T}$$ and $$\varepsilon > 0$$; when $$\varepsilon=0$$, the boundary conditions are applied only on the inflow part of $$\partial{\it\Omega}$$. For $$\varepsilon > 0$$ boundary layers are formed along $$x=1$$ and $$y=1$$. We compare the error estimator under an adaptive refinement algorithm, for $$\varepsilon = 1e-2, 1e-4, 0$$, where a bulk marking strategy is applied; in all cases, an initial uniform grid comprising $$512$$ simplices is used. The effectivities $$\eta/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_1$$ are shown in Fig. 7, where we notice that in all cases they remain bounded between $$1$$ and $$4$$: for $$\varepsilon = 1e-2$$, the effectivities are slowly increasing but would appear to be converging to an upper bound below $$4$$, while for $$\varepsilon = 0$$ the effectivities seem to converge to $$1$$. The results are remarkably consistent with those from Example 2. In Fig. 8, adaptively refined grids are shown for both $$\varepsilon = 1e-2$$ and $$\varepsilon = 1e-4$$. We notice that for $$\varepsilon = 1e-2$$ refinement has been carried out in the vicinity of the boundary layers but also along a line parallel to the convection $$\bf b$$; for $$\varepsilon = 1e-4$$, where the boundary layers are narrower, refinement has been carried out only at the boundaries. Fig. 7. View largeDownload slide Example 4: (a) Error bound effectivities for varying $$\varepsilon$$ Fig. 7. View largeDownload slide Example 4: (a) Error bound effectivities for varying $$\varepsilon$$ Fig. 8. View largeDownload slide Example 4: Adaptive mesh for $$\varepsilon = 1e-2$$, iteration number 13 (left) and for $$\varepsilon = 1e-4$$, iteration number 20 (right). Fig. 8. View largeDownload slide Example 4: Adaptive mesh for $$\varepsilon = 1e-2$$, iteration number 13 (left) and for $$\varepsilon = 1e-4$$, iteration number 20 (right). 6.6 Example 5 In our final numerical experiment, we show that the estimator can drive an adaptive strategy to resolve solutions with both interior and boundary conditions. To this end, we consider the same test problem as in Schötzau & Zhu (2009, Example 3). In this case, the problem is to find $$u$$ such that \begin{alignat*}{2} -\varepsilon{\it\Delta} u+{\bf b}\cdot \nabla u &= 0, &(x,y) \in {\it \Omega} = (-1,1)^2, \\ u &= 0 &\text{on } x=-1 \text{ and } y = 1, \\ u &= \tanh\left(\frac{1-y}{\epsilon}\right) &\text{on $x=1$}, \\ u &= \frac{1}{2}\left(\tanh\left(\frac{x}{\epsilon}\right)+1\right) &\text{on $y=-1$}, \end{alignat*} where $${\bf b} = (-\sin\frac{\pi}{6},\cos\frac{\pi}{6})^{\rm T}$$. The boundary conditions cause $$u$$ to have an internal layer along the line $$y+\sqrt{3}x = -1$$ and boundary layers at the outflow boundaries. Since there is no exact solution available for this problem, reference solutions for varying $$\varepsilon$$ have been obtained using a different a posteriori estimator (namely, the $$\eta_{\rm SZ}$$ estimator) within an adaptive framework and setting $$p=4$$ to ensure that the reference solution is of sufficient accuracy. In this case, we apply a fixed fraction strategy with $$25\%$$ of the elements chosen for refinement and $$10\%$$ chosen for derefinement at each iteration. For both $$\varepsilon = 10^{-2}$$ and $$\varepsilon = 10^{-4}$$, we begin with a uniform mesh comprising $$512$$ right-angled triangles, which are not aligned with the direction of flow, and set $$p=1$$. Figure 9 reports the error convergence history for both these cases, while Fig. 10 shows the resultant meshes after $$14$$ iterations of the adaptive algorithm for $$\varepsilon = 10^{-2}$$ and $$\varepsilon = 10^{-4}$$, respectively. Fig. 9. View largeDownload slide Example 5: Error convergence for adaptive strategy for $$\varepsilon=1e-2$$ and $$\varepsilon=1e-4$$. Fig. 9. View largeDownload slide Example 5: Error convergence for adaptive strategy for $$\varepsilon=1e-2$$ and $$\varepsilon=1e-4$$. Fig. 10. View largeDownload slide Example 5: Adaptive mesh for $$\varepsilon = 1e-2$$ (left) and for $$\varepsilon = 1e-4$$ (right) at iteration number $$14$$. Fig. 10. View largeDownload slide Example 5: Adaptive mesh for $$\varepsilon = 1e-2$$ (left) and for $$\varepsilon = 1e-4$$ (right) at iteration number $$14$$. The error indicator leads to refinement in the regions of both the boundary and interior layers for both values of $$\varepsilon$$. For $$\varepsilon=10^{-4}$$, a pre-asymptotic regime is observed, as expected; when the mesh is fine enough to resolve the boundary and interior layers, the error begins to converge at the same rate as for $$\varepsilon = 10^{-2}$$. References Ainsworth M. , Allendes A. , Barrenechea G. R. & Rankin R. ( 2013 ) Fully computable a posteriori error bounds for stabilized FEM approximations of convection-reaction-diffusion problems in three dimensions. Internat. J. Numer. Methods Fluids 73 , 765 – 790 . Arnold D. N. ( 1982 ) An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal. , 19 , 742 – 760 . Google Scholar CrossRef Search ADS Cohen A. , Dahmen W. & Welper G. ( 2012 ) Adaptivity and variational stabilization for convection-diffusion equations. ESAIM Math. Model. Numer. Anal. , 46 , 1247 – 1273 . Google Scholar CrossRef Search ADS Creusé E. & Nicaise S. ( 2013 ) A posteriori error estimator based on gradient recovery by averaging for convection-diffusion-reaction problems approximated by discontinuous Galerkin methods. IMA J. Numer. Anal. , 33 , 212 – 241 . Google Scholar CrossRef Search ADS Dahmen W. , Huang C. , Schwab C. & Welper G. ( 2012 ) Adaptive Petrov-Galerkin methods for first order transport equations. SIAM J. Numer. Anal. , 50 , 2420 – 2445 . Google Scholar CrossRef Search ADS Ern A. , Stephansen A. F. & Vohralík M. ( 2010 ) Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems. J. Comput. Appl. Math. , 234 , 114 – 130 . Google Scholar CrossRef Search ADS Georgoulis E. H. , Hall E. & Makridakis C. ( 2014 ) Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems . Barrett Lectures on Discontinuous Galerkin Methods. IMA Vol. Math. Appl. , vol. 157 , Cham : Springer , 195 – 207 . Georgoulis E. H. & Lasis A. ( 2006 ) A note on the design of $$hp$$-version interior penalty discontinuous Galerkin finite element methods for degenerate problems. IMA J. Numer. Anal. , 26 , 381 – 390 . Google Scholar CrossRef Search ADS Houston P. , Schwab C. & Süli E. ( 2002 ) Discontinuous $$hp$$-finite element methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. , 39 , 2133 – 2163 (electronic). Google Scholar CrossRef Search ADS Johnson C. & Pitkäranta J. ( 1986 ) An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math. Comput. , 46 , 1 – 26 . Google Scholar CrossRef Search ADS Makridakis C. & Nochetto R. H. ( 2006 ) A posteriori error analysis for higher order dissipative methods for evolution problems. Numer. Math. , 104 , 489 – 514 . Google Scholar CrossRef Search ADS Roos H.-G. , Stynes M. & Tobiska L. ( 2008 ) Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems , vol. 24 , 2nd edn . Springer Series in Computational Mathematics , Berlin : Springer . Sangalli G. ( 2001 ) A robust a posteriori estimator for the residual-free bubbles method applied to advection-diffusion problems. Numer. Math. , 89 , 379 – 399 . Google Scholar CrossRef Search ADS Sangalli G. ( 2008 ) Robust a-posteriori estimator for advection-diffusion-reaction problems. Math. Comp. , 77 , 41 – 70 (electronic). Google Scholar CrossRef Search ADS Schötzau D. & Wihler T. P. ( 2010 ) A posteriori error estimation for $$hp$$-version time-stepping methods for parabolic partial differential equations. Numer. Math. , 115 , 475 – 509 . Google Scholar CrossRef Search ADS Schötzau D. & Zhu L. ( 2009 ) A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations. Appl. Numer. Math. , 59 , 2236 – 2255 . Google Scholar CrossRef Search ADS Verfürth R. ( 2005 ) Robust a posteriori error estimates for stationary convection-diffusion equations. SIAM J. Numer. Anal. , 43 , 1766 – 1782 (electronic). Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

An a posteriori error bound for discontinuous Galerkin approximations of convection–diffusion problems

Loading next page...
 
/lp/ou_press/an-a-posteriori-error-bound-for-discontinuous-galerkin-approximations-T8Qp805f86
Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx065
Publisher site
See Article on Publisher Site

Abstract

Abstract An a posteriori bound for the error measured in the discontinuous energy norm for a discontinuous Galerkin (dG) discretization of a linear one-dimensional stationary convection–diffusion–reaction problem with essential boundary conditions is presented. The proof is based on a conforming recovery operator inspired from a posteriori error bounds for the dG method for first-order hyperbolic problems. As such, the bound remains valid in the singular limit of vanishing diffusion. Detailed numerical experiments demonstrate the independence of the quality of the a posteriori bound with respect to the Péclet number in the standard dG-energy norm, as well as with respect to the viscosity parameter. 1. Introduction The interaction between convection and diffusion in many physical systems often gives rise to multiscale solution behaviour, when the convection is the dominant phenomenon. These are typically manifested as boundary and/or interior layers or even as discontinuities, in cases of diffusion-degenerate problems. The stable discretization for convection–diffusion problems in the contexts of finite difference and finite element methods is now relatively well understood. However, such lower-dimensional solution features require carefully selected variable mesh resolution across the computational domain to be efficiently resolved. The location of such multiscale features may not be a priori available. This motivates the need for adaptive algorithms. In the context of finite element methods, there exists a growing literature on the derivation of adaptive algorithms for (steady-state and transient) convection–diffusion problems, based on local error estimation, either in an ad hoc fashion or via reliable a posteriori bounds. A key question in the respective literature is the independence (robustness in the terminology of most references) of the derived a posteriori estimators with respect to the, so-called, Péclet number, i.e., the ratio between the convection and the diffusion magnitudes. To achieve such error control with upper and lower a posteriori bounds, several approaches have been considered, with the most popular being the augmentation of the energy norm by a negative/fractional Sobolev norm of the skew-symmetric part of the differential operator (Verfürth, 2005; Sangalli, 2008; Schötzau & Zhu, 2009; Ern et al., 2010; Creusé & Nicaise, 2013); we refer to Sangalli (2008) for an insightful exposition of the differential operator’s stability properties with respect to the Péclet number. We also mention the use of subgrid problems (Sangalli, 2001), the recent adaptive variational stabilization methods (Cohen et al., 2012; Dahmen et al., 2012; Ainsworth et al., 2013), focusing on fully computable a posteriori bounds for the standard energy norm. Generally speaking, the aforementioned works (with the exception of Cohen et al., 2012; Dahmen et al., 2012, which introduce new stabilized numerical methods aiming at sharp error control) are based on adapting the a posteriori error analysis for the purely elliptic problem to the convection–diffusion(–reaction) problem yielding error estimates controlling the error with upper and lower bounds. As a result, although the a posteriori error bounds in Verfürth (2005), Sangalli (2008), Ern et al. (2010), Schötzau & Zhu (2009) and Creusé & Nicaise (2013) hold for each fixed positive diffusion parameter, the derived estimators are not applicable when the diffusion vanishes; this limits the possible applicability of such bounds to realistic nonlinear shock computations. Given the very successful computational performance of standard stabilized methods, such as discontinuous Galerkin (dG) methods with various fluxes, we would like to seek an alternative approach to their error control. Our work is motivated by the following issues: (1) whether it is possible to prove a posteriori bounds for singularly perturbed problems of this kind starting from a posteriori error control of the hyperbolic limit problem and (2) whether such a posteriori bounds for errors measured in standard energy-type norms (i.e., norms with respect to which the method’s bilinear forms are coercive) would hold with error bound constants that remain uniformly bounded with respect to vanishing viscosity and the Péclet number. Our answer to both questions is that, under certain conditions, such error bounds are possible. We stress that, in the present context of convection–diffusion problems, we address these questions with regard to the proof of upper bounds and their respective validity at the zero viscosity limit. In addition, we address, to a certain extent, similar questions related to lower bounds and we thoroughly investigate computationally the behaviour of the estimators. We derive a new a posteriori error bound for dG discretizations of a linear stationary convection–diffusion–reaction problem, utilizing ideas from error control for first-order hyperbolic problems. In fact, we show that the answers to the questions (1) and (2) above are positive, provided we have at our disposal a technique that leads to error control of the limiting hyperbolic problem. This can be done to the full generality considering the one-dimensional problem by utilizing the reconstructions in Makridakis & Nochetto (2006). The multidimensional hyperbolic problem was partially addressed in Georgoulis et al. (2014), and its full generality is an important open problem. Thus, we have chosen to present our results in full detail for the one-dimensional problem. It is clear from the analysis that similar results can be obtained in higher dimensions using the ideas in Georgoulis et al. (2014). Since such results are still not available in their full generality, we have chosen to consider numerical experiments to highlight the applicability of the proposed approach in higher dimensions. To this end, we derive a posteriori upper bounds for the dG method for the convection–diffusion–reaction problem, which remain valid in the hyperbolic limit of vanishing diffusion and at the same time remain bounded uniformly with respect to the Péclet number (see Theorems 4.2 and 4.3 and Section 4.2). The upper bounds are complemented with standard lower bounds of the flux jump term whose constant, as expected, degenerates with the Péclet number tending to infinity. To further theoretically justify the computationally observed robustness of the proposed a posteriori estimators with respect to the Péclet number, additional nonstandard lower bounds are also derived using classical asymptotic analysis results (cf. Theorem 5.2). The asymptotically optimal behaviour of the a posteriori error estimator, as well as its superiority on controlling the natural dG-energy norm of the problem, is observed in practice on a number of one-dimensional numerical examples. Moreover, we observe numerically that the quality of the derived a posteriori bounds with respect to the energy norm does not deteriorate with increasing Péclet number (when the viscosity tends to zero), as opposed to the respective estimators from Verfürth (2005), Sangalli (2008), Ern et al. (2010), Schötzau & Zhu (2009) and Creusé & Nicaise (2013). To this end, numerical comparisons with known estimators from the literature are performed. Finally, we consider some two-dimensional numerical experiments to highlight the performance of the estimators in higher dimensions. The remainder of this work is structured as follows. In Section 2, we present the model problem and the dG method for the one-dimensional case. In Section 3, we define and discuss the properties of the reconstruction used in the a posteriori analysis in Section 4. Section 4 contains the main results of this work, providing a posteriori bounds for various regimes of equation coefficients, while in Section 5 lower bounds are discussed. Section 6 contains a number of numerical experiments highlighting, in particular, the viscosity independence of the a posteriori error bounds with respect to the Péclet number, along with use of new a posteriori bounds within adaptive algorithms. 2. Preliminaries 2.1 Model problem and notation We shall make use of standard notation for the $$L^2$$-inner product $$(\cdot,\cdot)_{\omega}$$ on an interval $$\omega\subset \mathbb{R}$$, along with the respective $$L^2$$-norm $$\|{\cdot}\|_{\omega}$$; the subscript $$\omega$$ will be suppressed when $$\omega \equiv I:=(\alpha,\beta)\subset\mathbb{R}$$, the computational interval. The standard Hilbertian Sobolev spaces $$H^1(\omega)$$ and $$H^1_0(\omega)$$ will also be used. In $$I$$, we consider the boundary value problem \begin{equation}\label{pde} -\varepsilon u ''+b u'+cu =f\quad\text{in } I,\qquad u(\alpha) = 0 = u(\beta), \end{equation} (2.1) for $$f\in L^2(I)$$, $$\varepsilon > 0$$, $$0\ne b\in \mathbb{R}$$ and $$c:\bar{I}\to\mathbb{R}$$ a continuous function; in weak form this reads, find $$u\in H^1_0(I)$$, such that \begin{equation}\label{pde_weak} \varepsilon (u ',v')+(b u',v)+(cu,v) =(f,v)\quad\text{for all } v\in H^1_0(I). \end{equation} (2.2) We define the inflow node of $$I$$ by $$\partial_- I:=\alpha$$ if $$b>0$$ and $$\partial_- I:=\beta$$ if $$b<0$$; the respective outflow node is defined by $$\partial_+ I:=\{\alpha,\beta\}\backslash \partial_- I$$. The singular limit case $$\varepsilon =0$$ (together with the respective removal of the outflow boundary condition) will also be discussed; the respective hyperbolic problem in weak form reads, find $$u\in H^1_-(I)$$, such that \begin{equation}\label{pde_weak_hyp} (b u',v)+(cu,v) =(\,f,v)\quad\text{for all } v\in H^1_-(I), \end{equation} (2.3) with $$H^1_-(I):=\{v\in H^1(I): v(\partial_- I)=0\}$$. We note that completely analogous a posteriori bounds to the ones given below can be shown for mildly more general coefficients in (2.1) than those assumed above. We shall refrain from doing so in the interest of simplicity of the presentation, settling for providing comments in this direction at various parts of this work. Let $$\mathcal{T}$$ be a subdivision of the domain $$I$$ into elements $$T_i:=(x_{i-1},x_i)$$ for $$i=1,\dots, N$$, $$x_0:=\alpha$$, $$x_{N}:=\beta$$; each element $$T_i$$ has exactly one inflow node $$\partial_-T_i\in \{x_{i-1},x_i\}$$, depending on the sign of $$b$$. Also let $$\partial_+T_i:=\{x_{i-1},x_i\}\backslash \partial_-T_i$$. For notational convenience, we shall also denote by $$\mathcal{N}_-$$ (resp. $$\mathcal{N}_+$$) the set of indices of all interior nodes $$i=1,\dots, N-1$$, together with the index of the inflow node $$\partial_-I$$ (resp. outflow node $$\partial_+I$$). Further, we consider the corresponding (discontinuous) elementwise polynomial space \[ \mathbb{V}_p:=\left\{v\in L^2(I): v|_{(x_{i-1},x_{i})}\in\mathcal{P}_p(x_{i-1},x_{i})\right\}\!, \] where $$\mathcal{P}_p(x_{i-1},x_{i})$$ is the space of polynomials of degree at most $$p\in\mathbb{N}$$. For $$v\in \mathbb{V}_p$$, $$i=0,\dots,N$$, we define the upwind jump $$\lfloor{v}\rfloor(x_i)$$ across the node $$x_i$$ by $$\lfloor{v}\rfloor(x_i):=\lim_{\delta\to 0} ( v(x_i+b\delta) -v(x_i-b\delta) )$$, $$\delta>0$$, adopting the conventions $$\lfloor{v}\rfloor(\alpha) = v(\alpha)$$, $$\lfloor{v}\rfloor(\beta) = v(\beta)$$, i.e., the values taken from outside the domain $$I$$ are set to be equal to the homogeneous boundary conditions. We also define the average and jump across $$x_i$$ by \[ \{ \kern -1.6mm \{{v}\} \kern -1.6mm \}\ (x_i) := (v(x_i^+) +v(x_i^-))/2\quad\text{ and} \quad [ \kern -.7mm [{v}] \kern -.7mm ](x_i) := v(x_i^-) -v(x_i^+), \] respectively, with $$v(x_i^+)$$ and $$v(x_i^-)$$ denoting the traces from the right and from the left, respectively; we also adopt the conventions $$\{ \kern -1.6mm \{{v}\} \kern -1.6mm \}(\alpha)= v(\alpha)$$, $$[ \kern -.7mm [{v}] \kern -.7mm ](\alpha)= -v(\alpha)$$ and $$\{ \kern -1.6mm \{{v}\} \kern -1.6mm \}(\beta)=[ \kern -.7mm [{v}] \kern -.7mm ](\beta) = v(\beta)$$. Using this notation, we can also define the mesh-size function $$h:\bar{I}\to \mathbb{R}_+$$, given by $$h|_{T_i} = h_i:=x_i-x_{i-1}$$, $$i=1,\dots,N$$, $$h(\alpha) = h_1$$, $$h(\beta) = h_N$$, and $$h(x_i) := \{ \kern -1.6mm \{{h}\} \kern -1.6mm \}$$ for $$i=1,\dots,N-1$$. 2.2 The dG method We consider the (interior penalty) dG method reading, find $$u_h\in \mathbb{V}_p$$, such that \begin{equation}\label{dg_weak} \varepsilon B_{\rm d}(u_h,v_h)+ B_{\rm c}(u_h,v_h) = (\,f,v_h) \end{equation} (2.4) for all $$v_h\in \mathbb{V}_p$$, with \[ \begin{aligned} B_{\rm d}(w,v) := & \sum_{i=1}^N (w',v')_{T_i}-\sum_{i=0}^{N} \left(\{ \kern -1.6mm \{{w'}\} \kern -1.6mm \}[ \kern -.7mm [{v}] \kern -.7mm ]+ \{ \kern -1.6mm \{{ v'}\} \kern -1.6mm \}[ \kern -.7mm [{w}] \kern -.7mm ]-\sigma[ \kern -.7mm [{w}] \kern -.7mm ][ \kern -.7mm [{v}] \kern -.7mm ]\right)(x_i)\\ \end{aligned} \] and \[ B_{\rm c}(w,v):= \sum_{i=1}^N\Big( \big(bw'+cw,v\big)_{T_i} + |b|(\lfloor{w}\rfloor v^+)(\partial_-T_i)\Big), \] with $$v_h^+(\partial_-T_i)$$ denoting the value of $$v_h$$ on the inflow node $$\partial_-T_i$$ from within $$T_i$$, and $$\sigma> 0$$ given by $$\sigma(x_i) = C_{\sigma}/h(x_i)$$, $$i=0,\dots,N$$, for some user-defined constant $$C_{\sigma}>1$$. The corresponding dG-energy norm associated to $$B_{\rm d}$$ is defined by \[ | \kern -.25mm \|{w}| \kern -.25mm \| := \left(\sum_{i=1}^N \|{w'}\|_{T_i}^2 +\sum_{i=0}^{N}(\sigma[ \kern -.7mm [{w}] \kern -.7mm ]^2)(x_i)\right)^{\frac{1}{2}}. \] We also define the dG-energy norm for the method by \begin{equation}\label{en_norm} | \kern -.25mm \|{w}| \kern -.25mm \|_q : = \left(\varepsilon| \kern -.25mm \|{w}| \kern -.25mm \|^2 + \|{qw}\|_{}^2 + \frac{|b|}{2} \sum_{i=0}^{N}[ \kern -.7mm [{w}] \kern -.7mm ]^2(x_i)\right)^{\frac{1}{2}}, \end{equation} (2.5) with the subscript indicating the weight in the $$L^2$$-norm component of the dG-energy norm. The bilinear form $$B_{\rm d}$$ can be shown to be coercive in this norm in $$H^1_0\times H^1_0$$, and also in $$\mathbb{V}_p\times \mathbb{V}_p$$, provided the penalty constant $$C_{\sigma}$$ is chosen sufficiently large; for a proof we refer, e.g., to Arnold (1982) and Houston et al. (2002). For the bilinear form $$B_{\rm c}$$, the identity \[ B_{\rm c}(w,w) = \|{\sqrt{c}w}\|_{}^2 + \frac{|b|}{2} \sum_{i=0}^{N}[ \kern -.7mm [{w}] \kern -.7mm ]^2(x_i) \] holds for $$c\ge0$$ (see, e.g., Johnson & Pitkäranta, 1986; Houston et al., 2002). Weaker assumptions on $$c$$ will be discussed in Section 4.3. We note that the theory presented below also remains valid for nonsymmetric and incomplete interior penalty dG methods; this has not been carried through explicitly to minimize the notational overhead. 3. Reconstruction We begin by defining a reconstruction of the approximate solution, which is closely related to the optimal-order time reconstruction for dG time-stepping schemes, presented in Makridakis & Nochetto (2006). This will be a crucial ingredient in the proof of the a posteriori bound below, enabling its validity up to, and including, the singular limit $$\varepsilon =0$$. Definition 3.1 (Reconstruction). For each $$i=1,\dots,N$$, we define the dG reconstruction $$\hat{u}_h\in \mathbb{V}_{p+1}$$ of $$u_h\in \mathbb{V}_p$$ on each $$T_i$$, by the relations \begin{equation}\label{space_rec_def_formula} \big( (b\hat{u}_h)', v_h\big)_{T_i} = \big( (b u_h)', v_h\big)_{T_i} + |b|(\lfloor{u_h}\rfloor v_h^+)(\partial_-T_i) \end{equation} (3.1) for all $$v_h\in \mathbb{V}_p$$, and $$\hat{u}_h(\partial_-T_i) = u_h^-(\partial_-T_i)$$, with $$u_h^-(\partial_-T_i)$$ denoting the value of $$u_h$$ on the inflow node $$\partial_-T_i$$ from outside $$T_i$$; when $$\partial_-T_i=\partial_-I$$, we set $$\hat{u}_h(\partial_-T_i) = 0$$, i.e., the inflow boundary value. The following lemma shows the well-posedness of this reconstruction. The proof is essentially given in Makridakis & Nochetto (2006, Lemma 2.1) in the context of dG time-stepping methods; it is included here for completeness of the presentation. Lemma 3.2. The dG reconstruction $$\hat{u}_h$$ is uniquely defined and we have \begin{equation}\label{space_rec} \big( (b \hat{u}_h)' + cu_h, v_h\big) = B_{\rm c}(u_h,v_h) \end{equation} (3.2) for all $$v_h\in \mathbb{V}^n$$. Moreover, $$\hat{u}_h$$ is continuous in $$I$$. Proof. Identity (3.2) is evident from Definition 3.1. To show the continuity of the reconstruction for all but the outflow element, we integrate (3.1) by parts and we use the property $$\hat{u}_h(\partial_-T_i) = u_h^-(\partial_-T_i)$$, to obtain \begin{equation}\label{reczero} \int_{x_{i-1}}^{x_{i}} b(\hat{u}_h-u_h) v'_h\,\mathrm{d} x =|b|\left(\left(\hat{u}_h-u_h^+\right) v_h^+\right)(\partial_+T_i). \end{equation} (3.3) Now, setting $$v_h$$ to be a constant on $$T_i$$, we deduce $$\hat{u}_h(\partial_+T_i)=u_h^+(\partial_+T_i)$$, with $$u_h^+(\partial_+T_i)$$ denoting the value of $$u_h$$ at the outflow node $$(\partial_+T_i)$$ taken from within $$T_i$$. This implies continuity at all interior nodes $$x_i$$, $$i=1,\dots,N-1$$ and that the left-hand side of (3.3) is equal to zero for all $$v_h'\in\mathcal{P}_{p-1}(x_{i-1},x_i)$$. Hence, $$\pi_{p-1}(\hat{u}_h-u_h)=0$$ on $$(x_{i-1},x_i)$$, with $$\pi_{q}$$ denoting the orthogonal $$L^2$$-projection onto $$\mathbb{V}_q$$. This, together with the exactness at the interval end points, implies the uniqueness of the reconstruction on $$(x_{i-1},x_i)$$. Finally, (3.2) follows by summing for all $$i=1,\dots,N$$. □ A crucial stability property of the above reconstruction is the content of the following result. Lemma 3.3. For each element $$T_i$$, $$i=1,\dots,N$$, we have \[ \left \|{(\hat{u}_h-u_h)^{(s)} }\right\|_{T_i} \le C h_i^{1/2-s}|[ \kern -.7mm [{u_h}] \kern -.7mm ](\partial_-T_i)|, \] $$s=0,1$$, for a $$C>0$$ constant, independent of $$h_i$$ and of $$u_h\in\mathbb{V}_p$$, but dependent on the elemental polynomial degree $$p$$. Proof. The proof for $$s=0$$ follows completely analogously to the one of Makridakis & Nochetto (2006, Lemma 2.2) and is, therefore, omitted. The proof for $$s=1$$ then follows using the standard inverse estimate $$\|{(\hat{u}_h-u_h)' }\|_{T_i} \le Ch_i^{-1} \|{\hat{u}_h-u_h }\|_{T_i}$$ along with the bound for $$s=0$$. □ The dependence of the constant on the polynomial degree $$p$$ in Lemma 3.3 has been exactly understood (Schötzau & Wihler, 2010) but we refrain from retaining this dependence below in the interest of simplicity of presentation. Moreover, it is possible to extend the above to nonconstant $$b$$, provided $$b$$ does not change sign inside an element; this will be considered in detail elsewhere. Lemma 3.4. Let $$\bar{q}:I\to \mathbb{R}$$ be the elementwise constant function with $$\bar{q}|_{T_i}:=\sup_{x\in T_i}q(x)$$ for some non-negative function $$q:\bar{I}\to\mathbb{R}$$. The following bound holds: \begin{equation}\label{uhatmu_estimate2} \begin{aligned} | \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|_q^2 \le &\, \sum_{i\in \mathcal{N}_-} \left( C\left(\varepsilon h_i^{-1}+\bar{q}^2h_i\right)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) +\frac{|b|}{2}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) \right) \end{aligned} \end{equation} (3.4) for some constant $$C>0$$, independent of $$\varepsilon$$, $$b$$, $$c$$, $$q$$, $$h$$, $$\hat{u}_h$$ and of $$u_h$$, but dependent on $$p$$. Proof. This follows by Lemma 3.3 and properties $$\hat{u}_h(\partial_+I)=u_h(\partial_+I)$$ and $$\hat{u}_h(\partial_-I)=0$$. □ Remark 3.5. We note that the second term on the right-hand side of (3.4) can be removed from both the left- and right-hand sides of (3.4), without affecting the validity of this bound. 4. A posteriori error bound We present an a posteriori error bound for $$\varepsilon>0$$ whose proof remains valid at the singular limit $$\varepsilon=0$$ also; the singular limit will be discussed in Section 4.2. In the interest of clarity of presentation, we shall first consider the simpler case of nonvanishing reaction coefficient $$c$$ in (2.1) (Section 4.1), before moving to the general case of small or vanishing reaction in Section 4.3. A comparison with known a posteriori bounds from the literature is given in Section 6.4. To this end, we begin by considering an auxiliary boundary value problem. For $$\varepsilon>0$$, we set $$\tilde{u}: = u + g$$, where $$g:I\to \mathbb{R}$$ is a linear function on $$I$$ with $$g(\partial_-I) = 0$$ and $$g(\partial_+I) = u_h(\partial_+I)$$. Then, $$\tilde{u}\in H^1(I)$$ satisfies the boundary value problem \begin{equation}\label{pde_mod} -\varepsilon \tilde{u} ''+b \tilde{u}'+c\tilde{u} =f+bg'+cg\quad\text{in } I,\qquad \tilde{u}(\partial_-I) = 0, \qquad \tilde{u}(\partial_+I) = u_h(\partial_+I). \end{equation} (4.1) For the singular limit case $$\varepsilon=0$$, we consider only the inflow boundary condition $$\tilde{u}(\partial_-I)=0$$, and we set $$\tilde{u}:= u$$. Upon defining $$\rho :=\tilde{u}-\hat{u}_h$$, we observe that $$\rho\in H^1_0(I)$$ when $$\varepsilon>0$$, and $$\rho\in H^1(I)$$ when $$\varepsilon=0$$. 4.1 The case $$\varepsilon>0$$ and $$c>0$$ Let $$\pi:L^2(I)\to \mathbb{V}_p$$ denote the orthogonal $$L^2$$-projection onto the finite element space. Then, we have the following error equation for the reconstruction. Lemma 4.1 (Error equation for the reconstruction). Set $$z:=bg'+cg$$ for brevity. For $$\varepsilon>0$$, we have \begin{equation}\label{err_eq} \begin{aligned} \varepsilon (\rho',v') + B_{\rm c}(\rho,v) = (\,f-\pi f+z,v) - \varepsilon (\hat{u}'_h,v') +\varepsilon B_{\rm d}(u_h,\pi v) +(\pi(cu_h)-c\hat{u}_h,v) \end{aligned} \end{equation} (4.2) for all $$v\in H^1_0(I)$$. For $$\varepsilon =0$$, (4.2) holds for all $$v\in H^1(I)$$. Proof. From the continuity of the reconstruction $$\hat{u}_h$$ on $$I$$, we have \begin{equation}\label{err_one} \begin{aligned} \varepsilon (\rho',v') + (b\rho'+c\rho,v) &= (\,f+z,v) - \varepsilon (\hat{u}'_h,v')- (b\hat{u}'_h+c\hat{u}_h,v) \\ &= (\,f+z,v) - \varepsilon (\hat{u}'_h,v') - (b\hat{u}'_h,\pi v) -(c\hat{u}_h,v) , \end{aligned} \end{equation} (4.3) using (4.1) and the fact that $$\hat{u}'_h\in \mathbb{V}_p$$, respectively. Using (3.2) and (2.4), the third term on the right-hand side of (4.3) can be written as \begin{equation}\label{err_three} (b\hat{u}'_h,\pi v) = B_{\rm c}(u_h,\pi v) -(cu_h,\pi v) =-\varepsilon B_{\rm d}(u_h,\pi v) + (\,f,\pi v) -(cu_h,\pi v). \end{equation} (4.4) Applying (4.4) on (4.3), along with standard properties of the orthogonal $$L^2$$-projection, results in the right-hand side being equal to \[ (\,f-\pi f+z,v) - \varepsilon \left(\hat{u}'_h,v'\right) +\varepsilon B_{\rm d}(u_h,\pi v) +(\pi(cu_h)-c\hat{u}_h,v), \] which already implies the result. □ For the proof of the a posteriori bound for $$\varepsilon>0$$, the following extension of the bilinear form $$B_{\rm d}$$ from $$\mathbb{V}_p\times \mathbb{V}_p$$ to $$(H^1(I)+\mathbb{V}_p)\times ( H^1(I)+\mathbb{V}_p)$$, given by \[ \begin{aligned} B_{\rm d}(w,v) := & \sum_{i=1}^N (w',v')_{T_i}-\sum_{i=0}^{N} \left(\{ \kern -1.6mm \{{(\pi w)'}\} \kern -1.6mm \}[ \kern -.7mm [{v}] \kern -.7mm ]+ \{ \kern -1.6mm \{{(\pi v)'}\} \kern -1.6mm \}[ \kern -.7mm [{w}] \kern -.7mm ]-\sigma[ \kern -.7mm [{w}] \kern -.7mm ][ \kern -.7mm [{v}] \kern -.7mm ]\right)(x_i),\\ \end{aligned} \] will be useful. The extended $$B_{\rm d}$$ is both coercive and continuous in $$(H^1(I)+\mathbb{V}_p)\times ( H^1(I)+\mathbb{V}_p)$$ with respect to the $$| \kern -.25mm \|{\cdot}| \kern -.25mm \|$$-norm for sufficiently large penalty constant $$C_{\sigma}$$ (see, e.g., Georgoulis & Lasis, 2006 for a proof). Recalling the definition of $$\bar{q}:I\to \mathbb{R}$$ for a function $$q$$ from Lemma 3.4 (applied to $$c$$ below) we have the first of the main results of this work. Theorem 4.2. For $$\varepsilon>0$$ and $$c>0$$, the following a posteriori bound holds: \begin{align}\label{apost_final} | \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\sqrt{c}}^2& \le C\Bigg(\big\|{c_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\big\|_{}^2\nonumber\\ &\qquad\quad{} + \sum_{i=1}^{N-1} \!\! \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) +\sum_{i\in \mathcal{N}_-}\big((\varepsilon\sigma+ \bar{c} h)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2\big)(x_i)\nonumber\\ &\qquad\quad{} +((\varepsilon\sigma+ c_{\rm out} |b|)u_h^2)(\partial_+I)\Bigg) +\frac{|b|}{2}\sum_{i\in \mathcal{N}_-} [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i), \end{align} (4.5) with $$c_{\rm osc} := \min\{c^{-1/2}, \varepsilon^{-1/2}h\}$$, $$c_{\rm out}:= \max\{\|{cb^{-1}}\|_{},\|{bc^{-1}}\|_{}\}$$ and a constant $$C$$, independent of $$f$$, $$\varepsilon$$, $$b$$, $$c$$, $$h$$ and of $$u_h$$, but dependent on $$p$$ and on the size of the computational domain $$\beta-\alpha$$. Proof. Upon setting $$v=\rho$$ on (4.2), the terms on the left-hand side of (4.2) give \begin{equation} \varepsilon (\rho',v') + (b\rho'+c\rho,v) = \varepsilon \|{\rho'}\|_{}^2 + \|{\sqrt{c}\rho}\|_{}^2, \end{equation} (4.6) while the terms on the right-hand side of (4.2) can be bounded as follows. Since for $$\varepsilon>0$$, we have $$\hat{u}_h\in H^1(I)$$, with $$\hat{u}_h(\partial_-I) =0$$, $$\hat{u}_h(\partial_+I) =u_h(\partial_+I)$$ and $$\rho\in H^1_0(I)$$, we can deduce \[ (\hat{u}'_h,\rho') = B_{\rm d}(\hat{u}_h,\rho)+{\rm sign}(b) ((\pi \rho)'u_h)(\partial_+I), \] which implies, upon straightforward manipulation, \begin{align}\label{apost_one} B_{\rm d}(u_h,\pi\rho) -(\hat{u}'_h,\rho') & = B_{\rm d}(u_h,\pi\rho-\rho) + B_{\rm d}(u_h-\hat{u}_h,\rho)\nonumber\\ &\quad -{\rm sign}(b) ((\pi \rho)'u_h)(\partial_+I). \end{align} (4.7) For the first term on the right-hand side of (4.7), we set $$\tilde{v} =\pi\rho-\rho$$ (noting that $$\pi \tilde{v} =0$$ by construction) and we perform integration by parts on the elemental integrals: \[ \begin{aligned} B_{\rm d}(u_h,\tilde{v}) = & \sum_{i=1}^N (u'_h,\tilde{v}')_{T_i}-\sum_{i=0}^{N} \left(\{ \kern -1.6mm \{{u'_h}\} \kern -1.6mm \}[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]-\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ][ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]\right)(x_i)\\ = & \sum_{i=1}^{N-1} \left( [ \kern -.7mm [{u'_h\tilde{v}}] \kern -.7mm ]- \{ \kern -1.6mm \{{u'_h}\} \kern -1.6mm \}[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]+\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ][ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]\right)(x_i)+\sigma u_h\tilde{v}(x_N)+\sigma u_h\tilde{v}(x_0) , \end{aligned} \] since $$(u''_h,\tilde{v})_{T_i} =0$$ for all $$i=1,\dots,N$$, by the orthogonality of $$\pi$$. Further, noting the formula $$[ \kern -.7mm [{u'_h\tilde{v}}] \kern -.7mm ] = [ \kern -.7mm [{u'_h}] \kern -.7mm ]\{ \kern -1.6mm \{{\tilde{v}}\} \kern -1.6mm \}+\{ \kern -1.6mm \{{u'_h}\} \kern -1.6mm \}[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]$$ for all interior nodes, we arrive at \begin{equation}\label{diff_one} B_{\rm d}(u_h,\tilde{v}) = \sum_{i=1}^{N-1} \left([ \kern -.7mm [{u'_h}] \kern -.7mm ]\{ \kern -1.6mm \{{\tilde{v}}\} \kern -1.6mm \}+\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ][ \kern -.7mm [{\tilde{v}}] \kern -.7mm ]\right)(x_i)+\sigma u_h\tilde{v}(x_N)+\sigma u_h\tilde{v}(x_0). \end{equation} (4.8) Using the approximation properties of $$\pi$$, we have \[ |[ \kern -.7mm [{\tilde{v}}] \kern -.7mm ](x_i) |\le C \|{ \sqrt{h} \rho'}\|_{(x_{i-1},x_{i+1})},\quad |\{ \kern -1.6mm \{{\tilde{v}}\} \kern -1.6mm \}(x_i) |\le C \|{ \sqrt{h} \rho'}\|_{(x_{i-1},x_{i+1})}, \] and $$|\tilde{v}(x_0) |\le C \|{ \sqrt{h} \rho'}\|_{T_1}$$, $$|\tilde{v}(x_N) |\le C \|{ \sqrt{h} \rho'}\|_{T_N}$$, for some (generic) constant $$C>0$$, independent of $$\rho$$ and of $$h$$. Hence, the left-hand side of (4.8) can be estimated by \begin{equation}\label{diff_oneone} \begin{aligned} |B_{\rm d}(u_h,\tilde{v})| \le C \left(\sum_{i=1}^{N-1} \left(h|[ \kern -.7mm [{u'_h}] \kern -.7mm ]|^2\right)(x_i) +\sum_{i=0}^{N} \left(\sigma |[ \kern -.7mm [{u_h}] \kern -.7mm ]|^2\right)(x_i)\right)^{1/2}\|{\rho'}\|_{} \\ \end{aligned} \end{equation} (4.9) respectively, for some (generic) constant $$C>0$$, independent of $$\rho$$ and of $$h$$. The second term on the right-hand side of (4.7) can be estimated using the continuity of the extended bilinear form: \begin{equation}\label{diff_two} |B_{\rm d}(u_h-\hat{u}_h,\rho)|\le C | \kern -.25mm \|{u_h-\hat{u}_h}| \kern -.25mm \| | \kern -.25mm \|{\rho}| \kern -.25mm \|. \end{equation} (4.10) For the last term on the right-hand side of (4.7), using a standard inverse estimate along with the stability of the local $$L^2$$-projection in the $$H^1$$-seminorm, we deduce \begin{equation}\label{diff_threenew} | ((\pi \rho)'u_h)(\partial_+I)| \le C \|{\rho'}\|_{}(h^{-1/2}|u_h|)(\partial_+I). \end{equation} (4.11) Also, from the orthogonality of the $$L^2$$-projection, we have \begin{equation}\label{diff_three} \begin{aligned} &|(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| \\ &\quad{} = |((\,f-cu_h)-\pi (\,f-cu_h),\rho-\chi) -(c(\hat{u}_h-u_h),\rho)| \\ &\quad{} \le \|{(\,f-cu_h)-\pi (\,f-cu_h)}\|_{}\|{\rho-\chi}\|_{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\|{\sqrt{c}\rho}\|_{} \end{aligned} \end{equation} (4.12) for any $$\chi\in\mathbb{V}_p$$. The choice $$\chi = 0$$ results in \begin{align}\label{diff_threeone} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le \big(\|{c^{-1/2}((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}\nonumber\\ &\quad{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\big) \|{\sqrt{c}\rho}\|_{}, \end{align} (4.13) whereas the choice $$\chi = \pi \rho$$ results in \begin{align}\label{diff_threetwo} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le C \big(\|{h((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}\|{\rho'}\|_{}\nonumber\\ &\quad{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\|{\sqrt{c}\rho}\|_{}\big). \end{align} (4.14) Setting $$c_{\rm osc} = \min\{c^{-1/2}, \varepsilon^{-1/2}h\}$$, the last two bounds can be combined into \begin{align}\label{diff_threeall} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le C\big(\|{c_{osc}((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}\nonumber\\ &\quad{} +\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}\big){| \kern -.25mm \|{\rho}| \kern -.25mm \|}_{\sqrt{c}}. \end{align} (4.15) Finally, recalling that $$z=bg'+cg$$, we also have \begin{equation}\label{diff_fournew} |(z,\rho)| \le \|{c^{-1/2}z}\|_{}\|{\sqrt{c}\rho}\|_{} \le Cc_{\rm out} \sqrt{|b|}|u_h|(\partial_+I)\|{\sqrt{c}\rho}\|_{}. \end{equation} (4.16) Hence, using (4.9), (4.10), (4.11), (4.15) and (4.16) on (4.2) with $$v=\rho$$, we arrive at \begin{align}\label{apost_prefinal} {| \kern -.25mm \|{\rho}| \kern -.25mm \|}_{\sqrt{c}}^2 = \varepsilon \|{\rho'}\|_{}^2+\|{\sqrt{c}\rho}\|_{}^2 & \le C\Bigg(\|{c_{osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2\nonumber\\ &\qquad{} +\varepsilon| \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|^2+\|{\sqrt{c}(\hat{u}_h-u_h)}\|_{}^2\nonumber\\ &\qquad{} +\sum_{i=1}^{N-1} (\varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2)(x_i) +\sum_{i=0}^{N} (\varepsilon\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\nonumber\\ &\qquad{} + c_{\rm out} |b|u_h^2(\partial_+I)\Bigg). \end{align} (4.17) Also, we have \begin{equation}\label{tzi} \begin{aligned} | \kern -.25mm \|{g}| \kern -.25mm \|_{\sqrt{c}}^2 = &\ \varepsilon\|{g'}\|_{}^2+ \|{\sqrt{c}g}\|_{}^2+ \left(\left(\varepsilon\sigma +\frac{|b|}{2}\right) u_h^2\right)(\partial_+I)\\ \le &\ C((\varepsilon\sigma +c_{\rm out}|b|) u_h^2)(\partial_+I). \end{aligned} \end{equation} (4.18) Finally, observing the identity $$u-u_h = -g+\rho+\hat{u}_h-u_h,$$ the triangle inequality implies \[ | \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\sqrt{c}}\le | \kern -.25mm \|{g}| \kern -.25mm \|_{\sqrt{c}} +| \kern -.25mm \|{\rho}| \kern -.25mm \|_{\sqrt{c}}+| \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|_{\sqrt{c}}. \] The result readily follows by combining (4.17), (4.18) and the reconstruction stability estimates from Lemma 3.4. □ Theorem 4.2 above assumes that the reaction coefficient $$c$$ is positive. Moreover, the relative size between the convection and the reaction coefficients determines the magnitude of the constant $$c_{\rm out}$$ and, therefore, may adversely affect the constant of the a posteriori bound. It is possible to remove these restrictive effects, by considering a variant of Theorem 4.2; this will be the content of Section 4.3. 4.2 The singular limit case $$\varepsilon=0$$ Crucially, the bound (4.5) remains valid at the singular limit $$\varepsilon =0$$. An inspection of the proof of Theorem 4.2, noting carefully that in this case we have $$g=z=0$$, gives \[ \begin{aligned} & \|{\sqrt{c}(u-u_h)}\|_{}^2 + \frac{|b|}{2} \sum_{i\in\mathcal{N}_-}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) + \frac{|b|}{2} (u-u_h)^2(\partial_+I)\\ &\quad{} \le C\left( \|{c^{-1/2} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 +\sum_{i\in\mathcal{N}_-} (\bar{c} h[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i) \right) + \frac{|b|}{2} \sum_{i\in\mathcal{N}_-}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i). \end{aligned} \] In this case, we also trivially have an a posteriori bound for the $$L^2$$-norm of the error \begin{equation}\label{singular_limit_apost} \|{\sqrt{c}(u-u_h)}\|_{}^2 \le \, C\left( \|{c^{-1/2} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 +\sum_{i\in\mathcal{N}_-} (\bar{c} h[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\right)\!; \end{equation} (4.19) cf. Makridakis & Nochetto (2006), for the respective result in the context of dG time stepping. 4.3 The case of small or vanishing reaction $$c$$ Now, we shall remove the assumption on the magnitude of the reaction coefficient by considering a variant of Theorem 4.2, which will be shown using a modified test function. The use of such testing in a priori analysis of dG methods is classical (see, e.g., Johnson & Pitkäranta, 1986). Theorem 4.3. Let $$\varepsilon>0$$, and assume that $$\varepsilon$$, $$b$$, $$c$$ are such that $$c+|b|-\varepsilon>0$$. Further, we set $$\gamma:I\to\mathbb{R}_+$$ with $$\gamma:=(c+|b|-\varepsilon)^{1/2}$$. Then, the following a posteriori bound holds: \begin{align}\label{apost_final_tw} | \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}^2 & \le C\Bigg(\|{\tilde{c}_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2\nonumber \\ &\qquad{} + \sum_{i=1}^{N-1} \!\! \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) +\sum_{i\in \mathcal{N}_-}\big((\varepsilon\sigma+ \bar{\gamma}^2 h)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2\big)(x_i) \nonumber\\ &\qquad{} +((\varepsilon\sigma+ c_{\rm out,2} |b|)u_h^2)(\partial_+I)\Bigg) +\frac{|b|}{2}\sum_{i\in \mathcal{N}_-} [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i), \end{align} (4.20) with $$\tilde{c}_{\rm osc} := \min\{\gamma^{-1}, \varepsilon^{-1/2}h\}$$, $$c\rm _{out,2}:= \max\{\|{\gamma^2 b^{-1}}\|_{},\|{b\gamma^{-2}}\|_{}\}$$ and a constant $$C$$, independent of $$f$$, $$\varepsilon$$, $$b$$, $$c$$, $$h$$ and of $$u_h$$, but dependent on $$p$$ and on $$\beta-\alpha$$. Moreover, in the singular limit $$\varepsilon =0$$, we have \begin{equation}\label{singular_limit_apost_two} \|{\gamma(u-u_h)}\|_{}^2 \le \, C\left( \|{\gamma^{-1} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 +\sum_{i\in\mathcal{N}_-} (\bar{\gamma}^2 h[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\right)\!, \end{equation} (4.21) which is valid also in the case of the pure advection problem, i.e., when $$c=0$$. Proof. We begin by defining the weight function $$\vartheta(x):= e^{-{\rm sign}(b)x}$$, $$x\in I$$ and we observe the trivial identity $$\vartheta'=-{\rm sign}(b)\vartheta$$. Further, we define the $$\vartheta$$-weighted $$L^2$$-norm by $$\|{(\cdot)}\|_{\vartheta}:=(\int_I\vartheta(x)(\cdot)^2(x)\,\mathrm{d} x)^{1/2}$$. Upon setting $$v=\rho \vartheta$$, on (4.2), and noting that $$v'= \rho'\vartheta + \rho\theta'=\rho'\vartheta-{\rm sign}(b)\rho\vartheta$$, the terms on the left-hand side of (4.2) give \begin{align} \varepsilon (\rho',v') + (b\rho'+c\rho,v) & = \varepsilon \|{\rho'}\|_{\vartheta}^2 + \|{\sqrt{c}\rho}\|_{\vartheta}^2 +(b-{\rm sign}(b) \varepsilon)(\rho',\rho\vartheta)\nonumber\\ & = \varepsilon \|{\rho'}\|_{\vartheta}^2 + \|{\gamma \rho}\|_{\vartheta}^2, \end{align} (4.22) upon integration by parts. The terms on the right-hand side of (4.2) can be bounded in the spirit of the proof of Theorem 4.2; we shall briefly repeat the estimates below focusing on the differences. Noting that $$v=\rho\vartheta\in H^1_0(I)$$, we have, as before, \begin{align}\label{apost_one_two} B_{\rm d}(u_h,\pi v) -(\hat{u}'_h,v') & = B_{\rm d}(u_h,\pi v-v) + B_{\rm d}(u_h-\hat{u}_h,v)\nonumber\\ &\quad{} -{\rm sign}(b) ((\pi v)'u_h)(\partial_+I). \end{align} (4.23) The first term on the right-hand side of (4.23), setting $$\tilde{v} =\pi v-v$$, performing integration by parts and using the approximation properties of $$\pi$$, as before, can be estimated by \begin{equation}\label{diff_oneone_two} |B_{\rm d}(u_h,\tilde{v})| \le C \left(\sum_{i=1}^{N-1} (h|[ \kern -.7mm [{u'_h}] \kern -.7mm ]|^2)(x_i) +\sum_{i=0}^{N} (\sigma |[ \kern -.7mm [{u_h}] \kern -.7mm ]|^2)(x_i)\right)^{1/2}\|{\rho'}\|_{\vartheta}, \end{equation} (4.24) observing the straightforward bound $$\|{v'}\|_{}\le C (\|{\rho'}\|_{\vartheta}+\|{\rho}\|_{\vartheta})\le C \|{\rho'}\|_{\vartheta}$$, with $$C$$ depending only on the interval $$I$$. The second term on the right-hand side of (4.23) can be estimated using the continuity of the extended bilinear form: \begin{equation}\label{diff_two_two} |B_{\rm d}(u_h-\hat{u}_h,v)|\le C | \kern -.25mm \|{u_h-\hat{u}_h}| \kern -.25mm \|\|{\rho'}\|_{\vartheta} , \end{equation} (4.25) since $$| \kern -.25mm \|{v}| \kern -.25mm \|_d=\|{v'}\|_{}\le C \|{\rho'}\|_{\vartheta}$$. For the last term on the right-hand side of (4.23), using a standard inverse estimate, we deduce \begin{equation}\label{diff_threenew_two} |((\pi v)'u_h)(\partial_+I)| \le C \|{v'}\|_{}(h^{-1/2}|u_h|)(\partial_+I)\le C \|{\rho'}\|_{\vartheta}(h^{-1/2}|u_h|)(\partial_+I). \end{equation} (4.26) Also, from the orthogonality of the $$L^2$$-projection, and working as before, we can have \begin{align}\label{diff_threeall_two} |(\,f-c\hat{u}_h-\pi (\,f-cu_h),\rho)| & \le C\big(\|{\tilde{c}_{\rm osc}((\,f-cu_h)-\pi (\,f-cu_h))}\|_{\vartheta}\nonumber\\ &\quad{} +\|{\gamma(\hat{u}_h-u_h)}\|_{\vartheta}\big)( \varepsilon \|{\rho'}\|_{\vartheta}^2+ \|{\gamma \rho}\|_{\vartheta}^2)^{1/2}, \end{align} (4.27) with $$\tilde{c}_{osc} = \min\{\gamma^{-1/2}, \varepsilon^{-1/2}h\}$$. Finally, we also have \begin{equation}\label{diff_fournew_two} |(z,v)| \le \|{\gamma^{-1}z}\|_{\vartheta}\|{\gamma \rho}\|_{\vartheta} \le Cc_{\rm out,2}\sqrt{|b|}|u_h|(\partial_+I)\|{\gamma\rho}\|_{\vartheta}. \end{equation} (4.28) Using the above bounds, we can arrive at \begin{align}\label{apost_prefinal_two} \varepsilon \|{\rho'}\|_{\vartheta}^2+\|{\gamma\rho}\|_{\vartheta}^2 & \le C\Bigg(\|{\tilde{c}_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2\nonumber\\ &\quad{} +\varepsilon| \kern -.25mm \|{\hat{u}_h-u_h}| \kern -.25mm \|^2+\|{\gamma(\hat{u}_h-u_h)}\|_{}^2\nonumber\\ &\quad{} +\sum_{i=1}^{N-1} (\varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2)(x_i) +\sum_{i=0}^{N} (\varepsilon\sigma[ \kern -.7mm [{u_h}] \kern -.7mm ]^2)(x_i)\nonumber\\ &\quad{} +c_{\rm out,2} |b|u_h^2(\partial_+I)\Bigg). \end{align} (4.29) Also, analogously to (4.18), we have $$| \kern -.25mm \|{g}| \kern -.25mm \|_{\gamma}^2 \le C((\varepsilon\sigma +c_{\rm out,2}|b|) u_h^2)(\partial_+I)$$. Finally, as before, the triangle inequality and the reconstruction stability estimates from Lemma 3.4 provide the first bound. The bound in the singular limit $$\varepsilon =0$$ follows by an inspection of the proof above with $$\varepsilon=g=z=0$$. □ Remark 4.4. The above result improves upon Theorem 4.2, in terms of the dependence with respect to the relative size between the reaction and the convection coefficients. Indeed, an $$L^2$$-component remains positive when $$c=0$$ and the constant $$c\rm _{out,2}$$ does not degenerate as $$c\to 0$$. 4.4 On extension to two space dimensions Here, we comment on the current obstacles on extending the above analysis to two (or more) space dimensions. First and foremost, the results for the respective hyperbolic problem presented in Georgoulis et al. (2014) are quite preliminary at this point and are valid only for special structured meshes in sufficient generality for the above analysis to carry through. A second challenge, of a rather technical nature, is the extension of the construction (4.1) in two dimensions. 5. On lower bounds The above a posteriori bound (4.5) can be complemented by the following lower bound, Theorem 5.1. This is derived by utilizing standard techniques and, as expected, does not provide enough information for very small values of $$\varepsilon$$ and moderate values of $$h.$$ A more detailed analysis using the structure of the possible boundary layer yields the second result in this section, Theorem 5.2, which establishes the quality of the estimator for small values of $$\varepsilon$$. Theorem 5.1. For each element $$T_i\in\mathcal{T}$$, $$i=1,\dots,N$$, we have \begin{align}\label{second_lower} \varepsilon^2 h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & \le C\left(h^2\tilde{{\rm Osc}}_i^2 + \varepsilon\|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + |b|^2 \|{u-u_h}\|_{T_i\cup T_{i+1}}^2 \right.\nonumber \\ &\quad{} \left.+|b|^2 \left( h_i [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_{i-1})+h_{i+1} [ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_{i})\right) + c_i \|{h\sqrt{c}(u-u_h)}\|_{T_i\cup T_{i+1}}^2\right)\!, \end{align} (5.1) with $${\rm Osc}_i: = \|{f-cu_h-\pi (\,f-cu_h)}\|_{T_i}$$, $$i=1,\dots,N$$ and $$\tilde{{\rm Osc}}_i :=({\rm Osc}_i^2+{\rm Osc}_{i+1}^2)^{1/2}$$, and $$C>0$$ a constant dependent only on $$p$$ and on the local quasi-uniformity of the mesh around $$x_i$$. Proof. We start by defining an elemental bubble function $$\beta_i$$ to be a quadratic polynomial on the element $$T_i$$, such that $$\beta_i(x_{i-1}) = 0 = \beta_i(x_i)$$ and $$\beta_i((x_{i-1}+x_i)/2) = 1$$, resulting in $$0\le \beta_i\le 1$$ and $$\|{\beta_i}\|_{T_i}^2 = C h_i$$, for a constant $$C$$ independent of $$h_i$$. For $$\psi : = h_i^2\phi\beta_i$$, with $$\phi:= \pi f +\varepsilon u_h'' -b\hat{u}_h'-\pi(c u_h)$$, we have, respectively, \begin{align} & \int_{T_i} \varepsilon (u-u_h)'\psi' -b(u-\hat{u}_h)\psi'+c(u-u_h)\psi \,\mathrm{d} x\nonumber\\ &\quad{} = \int_{T_i}\phi\psi \,\mathrm{d} x - \int_{T_i}(\,f-cu_h-\pi(\,f-cu_h))\psi \,\mathrm{d} x, \end{align} (5.2) which gives \begin{align} h_i^2\|{\sqrt{\beta_i} \phi}\|_{T_i}^2 & \le h_i^2 {\rm Osc}_i\|{\phi}\|_{T_i} + C \varepsilon h_i \|{(u-u_h)'}\|_{T_i}\|{\phi}\|_{T_i}\nonumber\\ &\quad{} + C |b| h_i\|{u-\hat{u}_h}\|_{T_i}\|{\phi}\|_{T_i}\nonumber\\ &\quad{} + \sqrt{c_i} h_i^2 \|{\sqrt{c}(u-u_h)}\|_{T_i}\|{\phi}\|_{T_i}, \end{align} (5.3) and, finally, using the (standard) bound $$\|{\phi}\|_{T_i} \le C \|{\sqrt{\beta_i} \phi}\|_{T_i}$$, for some $$C$$ independent of $$\phi$$ and of $$T_i$$, we get \begin{equation}\label{first_lower} \begin{aligned} h_i^2\|{ \phi}\|_{T_i}^2 \le &\ C\left( h_i^2 {\rm Osc}_i^2+ \varepsilon \|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i}^2 \right.\\ &\left.+ |b|^2\|{u-\hat{u}_h}\|_{T_i}^2 + c_i h_i^2 \|{\sqrt{c}(u-u_h)}\|_{T_i}^2\right)\!. \end{aligned} \end{equation} (5.4) Having estimated $$\|{ \phi}\|_{T_i}$$ by a norm of the error, we shall now estimate the term $$\varepsilon[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2$$ from above by the same norm of the error. To this end, we define the bubble function $$\tilde{\beta}_i$$ to be a linear polynomial on each $$T_i$$ and $$T_{i+1}$$, such that $$\tilde{\beta}_i(x_{i-1}) = 0 = \tilde{\beta}_i(x_{i+1})$$ and $$\tilde{\beta}_i(x_i) = 1$$ on $$T_i$$, resulting in $$0\le \tilde{\beta}_i\le 1$$ and $$\|{\tilde{\beta}_i}\|_{T_j}^2 = h_j/2$$, $$j=i,i+1$$. Now, for $$\tilde{\psi} : = \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\tilde{\beta}_i$$, we have, respectively, \begin{align} & \int_{T_i\cup T_{i+1}} \!\!\varepsilon (u-u_h)'\tilde{\psi}' -b(u-\hat{u}_h)\tilde{\psi}'+c(u-u_h)\tilde{\psi} \,\mathrm{d} x\nonumber \\ &\quad{} = \int_{T_i\cup T_{i+1}}\!\!\! \phi\tilde{\psi} \,\mathrm{d} x - \varepsilon([ \kern -.7mm [{u'_h}] \kern -.7mm ]\tilde{\psi})(x_i), \end{align} (5.5) noting that $$\int_{T_i\cup T_{i+1}}\!(\,f-cu_h-\pi(\,f-cu_h))\tilde{\psi} \,\mathrm{d} x =0$$, because $$\tilde{\psi} \in \mathbb{V}_p$$. This gives \begin{align} \varepsilon^2 h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & \le \|{\phi}\|_{T_i\cup T_{i+1}}\varepsilon h^{3/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\nonumber\\ &\quad{} + C\|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i\cup T_{i+1}}\varepsilon^{3/2} h^{1/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\nonumber\\ &\quad{} + C |b| \|{u-\hat{u}_h}\|_{T_i\cup T_{i+1}}\varepsilon h^{1/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i)\nonumber\\ &\quad{} + \sqrt{c_i} \|{\sqrt{c}(u-u_h)}\|_{T_i\cup T_{i+1}}\varepsilon h^{3/2}[ \kern -.7mm [{u'_h}] \kern -.7mm ](x_i), \end{align} (5.6) and, thus, \begin{align}\label{second_lower2} \varepsilon^2 h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & \le C\left( \|{h\phi}\|_{T_i\cup T_{i+1}}^2 + \varepsilon\|{\sqrt{\varepsilon}(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 \right.\nonumber\\ &\quad{} \left.+ |b|^2 \|{u-\hat{u}_h}\|_{T_i\cup T_{i+1}}^2 + c_i \|{h\sqrt{c}(u-u_h)}\|_{T_i\cup T_{i+1}}^2\right)\!. \end{align} (5.7) Using (5.4) to estimate the first term on the right-hand side of (5.7), and using Lemma 3.3, we finally arrive at the result. □ We remark that the remaining terms in the a posteriori bound above either involve jump-type terms of the numerical solution and are present also in the dG-energy norm or represent data oscillation. Indeed, it is possible to add the remaining jump-type terms from Theorems 4.2 and 4.3 to both sides of (5.1) to arrive at a complete lower bound. We note that the constant in the above lower bound is proportional to $$\varepsilon$$. This is to be expected as the continuity constant of the bilinear form $$\epsilon B_{\rm d} +B_{\rm c}$$ with respect to the dG-energy norm $$| \kern -.25mm \|{\cdot}| \kern -.25mm \|_{\gamma}$$ depends unfavourably on the Péclet number. On the other hand, as we shall see in the numerical experiments below, the term $$\varepsilon h [ \kern -.7mm [{u_h'}] \kern -.7mm ]^2$$ does not appear to be dominant in the pre-asymptotic regime (i.e., when $$h> \varepsilon$$) at least for problems involving boundary and interior layers. This observation can be further substantiated as follows. Let $$u_0\in H^1(I)$$ be the exact solution to the reduced problem $$bu'_0+cu_0 = f$$, $$u_0(\partial_- I)=0$$, and define $$s_0\in C^1(\bar{I})$$ to be a $$C^1$$-piecewise polynomial interpolant of $$u_0$$ (e.g., Hermite interpolant) with nodes $$x_i$$, $$i=0,1\dots,N$$. Also, let $$\pi_0:L^2(I)\to \mathbb{V}_0$$ denote the orthogonal $$L^2$$-projection operator onto constants on each $$T_i$$. We then have \begin{align}\label{lower_bound_two} \varepsilon h_i[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) & =\varepsilon h_i[ \kern -.7mm [{(s_0-u_h)'}] \kern -.7mm ]^2(x_i)\nonumber\\ &\le C\varepsilon \sum_{j=i,i+1}\|{(s_0-u_h)'}\|_{T_j}^2\nonumber\\ &=C\varepsilon \sum_{j=i,i+1}\|{(s_0-u_h-\pi_0(u-u_h))'}\|_{T_j}^2\nonumber\\ &\le C\varepsilon \sum_{j=i,i+1}h_j^{-2}\|{s_0-u_h-\pi_0(u-u_h)}\|_{T_j}^2\nonumber\\ &\le C\varepsilon \sum_{j=i,i+1}h_j^{-2}\Big( \|{u-u_h-\pi_0(u-u_h)}\|_{T_j}^2+ \|{u-u_0}\|_{T_j}^2+ \|{u_0-s_0}\|_{T_j}^2\Big)\nonumber\\ &\le C \sum_{j=i,i+1}\Big( \varepsilon \|{(u-u_h)'}\|_{T_j}^2+ \frac{\varepsilon }{h_j}\|{u-u_0}\|_{L^{\infty}(T_j)}^2+ \frac{\varepsilon }{h_j^{2}}\|{u_0-s_0}\|_{T_j}^2\Big), \end{align} (5.8) using inverse estimates and the local quasi-uniformity of the mesh, respectively. Now, classical asymptotic expansions (see, e.g., Roos et al., 2008, Theorem I.1.4) ensure that there exists a constant $$C>0$$, independent of $$\varepsilon$$, such that \begin{equation}\label{asymptotic} \|{u-u_0}\|_{L^{\infty}(I_{\zeta})}\le C\varepsilon, \end{equation} (5.9) with $$I_{\zeta}:=\{x\in I: {\rm dist} (x, \partial_+I)>\zeta\}$$, for fixed $$\zeta>0$$, sufficiently large, so that it does not include the boundary layer region of size $$\mathcal{O}(\varepsilon)$$ in the vicinity of $$\partial_+I$$. Also, since $$\|u\|_{L^{\infty}(I)}\le C$$ uniformly with respect to $$\varepsilon$$ (see, e.g., Roos et al., 2008, Theorem I.1.4), we have \[ \|{u-u_0}\|_{T_j\cap (I\backslash I_\zeta)}^2\le C|T_j\cap (I\backslash I_\zeta)|\le Ch_j \] for all $$T_j$$ with $$T_j\cap (I\backslash I_\zeta)\ne \emptyset$$, with $$C>0$$ a constant independent of $$\varepsilon$$ and of $$h_j$$. (Note that taking $$\varepsilon$$ small enough, we can have only one element $$T_j$$ such that $$T_j\cap (I\backslash I_\zeta)\ne \emptyset$$ in the mesh.) Moreover, since $$u_0\in H^1(I)$$, with $$u_0(x_i)=s_0(x_i)$$ from the interpolation property of $$s_0$$ at the nodes, the Poincaré–Friedrichs inequality implies \begin{equation}\label{pf_reduced} \|{u_0-s_0}\|_{T_i}\le \frac{h_i}{\pi} \|{(u_0-s_0)'}\|_{T_i} \end{equation} (5.10) for $$i=1,\dots,N$$. Furthermore, selecting $$s_0$$ to be exactly the Hermite interpolant of $$u_0$$ at the nodes $$x_i$$ and, assuming also $$f,c\in H^1(T_i)$$, we have $$u_0\in H^2(T_i)$$ and thus, \begin{equation}\label{pf_reduced_two} \|{u_0-s_0}\|_{T_i}\le \frac{h_i^2}{\pi^2} \|{(u_0-s_0)''}\|_{T_i}\le C h_i^2, \end{equation} (5.11) for some constant $$C>0$$ (which is, by construction, independent of $$\varepsilon$$). Combining the above estimates, we have the following result. Theorem 5.2. Let $$u_0\in H^{1+r}(I)$$, for $$r= 0,1$$. Then, for each element $$T_i\in\mathcal{T}$$, $$i=1,\dots,N$$ with $$T_i\in I_\zeta$$, we have \begin{equation}\label{lower_two_thm} \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon^3 h_i^{-1} + \varepsilon h_i^{2r} \big), \end{equation} (5.12) with the constant $$C>0$$, independent of $$\varepsilon$$, $$h$$ and of $$u$$, but possibly dependent on $$b$$, $$c$$ and on $$f$$. Moreover, we have \begin{equation}\label{lower_two_thm_two} \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon h_i^{-1} + \varepsilon h_i^{2r} \big) \end{equation} (5.13) for $$(T_i\cup T_{i+1}) \cap (I\backslash I_\zeta)\ne \emptyset$$, also. □ Now, assume for the moment that $$h_i>\varepsilon$$ for all $$i=1,\dots,N$$, and let $$k$$ be the index of the element $$T_k$$ having the outflow boundary $$\partial_+I$$ as one of its end points (i.e., $$k\in\{1,N\}$$) and that $$\varepsilon$$ is small enough, so that choosing $$\zeta =h_k$$ ensures the validity of (5.9) for this $$I_\zeta$$. In this setting, Theorem 5.2 implies \[ \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon h_i^r \big) \] for $$T_i\cup T_{i+1} \subset I_\zeta$$, and \begin{equation}\label{boundary_layer_element} \varepsilon h [ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \le C \big( \varepsilon\|{(u-u_h)'}\|_{T_i\cup T_{i+1}}^2 + \varepsilon h_i^{-1} \big) \end{equation} (5.14) for $$(T_i\cup T_{i+1}) \cap (I\backslash I_\zeta)\ne \emptyset$$, with $$k\in \{i,i+1\}$$. Therefore, fixing $$h$$ and letting $$\varepsilon\to 0$$ highlights the robustness of the proposed a posteriori estimator in the pre-asymptotic regime $$h_k>\varepsilon$$ for all elements in $$I_\zeta$$. Moreover, for $$\varepsilon$$ small enough, we can safely assume that $$T_k=I\backslash I_\zeta$$. For this element, the lower bound (5.14) shows that the estimator remains bounded as $$h_k\to \varepsilon$$. On the other hand, the good quality of the proposed estimator follows from Theorem 5.1 in the asymptotic regime $$h_k\le \varepsilon$$. In this latter case, dividing both sides of (5.1) by $$\varepsilon$$ and noticing that the last three terms on the right-hand side of (5.1) are of higher order with respect to the local mesh size, the boundedness of the local mesh-Péclet number $$|b|h_j/\varepsilon$$, $$j=i,i+1$$ ensures the good behaviour of the constant in (5.1). 6. Numerical experiments We present a number of examples to verify the quality of the estimators proposed in Theorems 4.2 and 4.3, especially their stability as the Péclet number $$|b|/\varepsilon\to \infty$$ and their ability to drive automatic adaptive mesh-refinement strategies. In Section 6.3, we discuss the a posteriori error indicator presented above, in conjunction with the augmented-norm-type estimators and in particular with the ones from Schötzau & Zhu (2009) and Ern et al. (2010). 6.1 Example 1 We begin by studying how the new a posteriori error estimate behaves for a problem with a smooth solution for various values of $$\varepsilon$$ under uniform mesh refinement. In this case, we choose $$(\alpha,\beta) = (0,1)$$, $$b=1$$, $$c=1$$ and pick $$f$$ such that $$u (x) = \sin(8\pi x)$$. Figure 1(a) shows the convergence of the error measured in the dG norm as the mesh is refined, while Fig. 1(b) shows the effectivities $$\eta/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ as the mesh is refined, for $$\varepsilon = 1,1e-1,\ldots,1e-4$$ and for $$\varepsilon = 0$$, with $$\eta$$ denoting the right-hand side of (4.20) without the constant, namely, \[ \begin{aligned} \eta:=&\|{\tilde{c}_{\rm osc} ((\,f-cu_h)-\pi (\,f-cu_h))}\|_{}^2 + \sum_{i=1}^{N-1} \!\! \varepsilon h[ \kern -.7mm [{u'_h}] \kern -.7mm ]^2(x_i) \\ &+\sum_{i\in \mathcal{N}_-}\left(\left(\varepsilon\sigma+ \bar{\gamma}^2 h+\frac{|b|}{2}\right)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2\right)(x_i) +((\varepsilon\sigma+ c_{\rm out,2} |b|)u_h^2)(\partial_+I). \end{aligned} \] Fig. 1. View largeDownload slide Example 1. Smooth problem. (a) Error convergence under uniform refinement. (b) Effectivity of estimator under uniform refinement. Fig. 1. View largeDownload slide Example 1. Smooth problem. (a) Error convergence under uniform refinement. (b) Effectivity of estimator under uniform refinement. As expected from a priori error analysis (see, e.g., Houston et al., 2002) and noting that, for a uniform mesh, $$h \sim N^{-1}$$, Fig. 1(a) shows that for $$\varepsilon > 0$$, convergence is $$\mathcal{O}(h)$$ in the asymptotic regime, whereas for $$\varepsilon = 0$$, the order is $$\mathcal{O}(h^{3/2})$$. We note that, for small $$\varepsilon$$, the jump terms in the dG-energy norm dominate on coarse meshes, hence initially $$\mathcal{O}(h^{3/2})$$ convergence is witnessed. For all $$\varepsilon$$, the effectivities remain bounded between $$0.9$$ and $$3.5$$, showing robustness of the error estimator as $$\varepsilon \rightarrow 0$$. In fact, as the mesh is refined, effectivities tend to a little over $$3$$ for $$\varepsilon >0$$, whereas, effectivities tend to $$1$$ for $$\varepsilon = 0$$. 6.2 Example 2 In this second example, we again select $$b=1$$ and $$c=1$$, but let $$f=1$$ on the domain $$(\alpha,\beta) = (0,1)$$; the exact solution exhibits a boundary layer near $$x=1$$, which narrows as $$\varepsilon$$ is reduced. Once again, we test the robustness of the estimator $$\eta$$, i.e., the right-hand side of (4.20), but this time we do so within an adaptive strategy. Elements are marked for refinement using a bulk criterion, where those elements contributing to 50% of the total error are refined; no elements are chosen for coarsening. In all cases, the computation starts from a uniform mesh comprising $$8$$ elements. We test with $$\varepsilon = 1, 1e-1,\ldots,1e-7$$ and $$\varepsilon = 0$$, noting there is no boundary layer in the latter case. Figure 2(a) shows the effectivities $$\eta/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ as the mesh is adaptively refined. In all cases, effectivities remain bounded between $$1$$ and $$3.5$$ and converge to a value just over $$3$$ for $$\varepsilon > 0$$ and to $$1$$ for $$\varepsilon = 0$$, highlighting the robustness of the proposed estimator. Fig. 2. View largeDownload slide Example 2. Boundary layer problem under adaptive refinement. (a) Effectivity of the error estimator. (b) Error convergence under uniform and adaptive refinement. Fig. 2. View largeDownload slide Example 2. Boundary layer problem under adaptive refinement. (a) Effectivity of the error estimator. (b) Error convergence under uniform and adaptive refinement. Figure 2(b) compares the error convergence for the adaptive strategy against uniform refinement for $$\varepsilon = 1e-2$$ and $$\varepsilon = 1e-5$$. For $$\varepsilon = 1e-2$$, on the latter meshes, the adaptive strategy gives well over an order of magnitude improvement in error over uniform refinement for a comparable number of degrees of freedom. For $$\varepsilon = 1e-5$$ and for around 5,000 degrees of freedom, the uniform strategy has failed to resolve the boundary layer successfully, while the adaptive strategy is approximately $$3$$ orders of magnitude superior for comparable degrees of freedom. 6.3 Example 3 We now test the bound from Theorem 4.3 by selecting $$c=0$$. As before, we let $$(\alpha,\beta) = (0,1)$$, $$b=1$$ but choose $$f = {\rm exp}(-10,000(x-0.5)^2)$$, which introduces an interior layer into the solution, narrowing as $$\varepsilon$$ is reduced. To prevent a boundary layer forming, we choose a homogeneous Neumann boundary condition on the right-hand side of the domain $$u^\prime(1) = 0$$; the estimator from Theorem 4.3 is modified in a trivial fashion to allow for this. Once again, we start with a uniform mesh of $$8$$ elements and carry out refinement using the same bulk criterion marking strategy as in Example 2, in this case for $$\varepsilon = 1e-1,\dots,1e-7$$. Figure 3(a) shows the effectivities as the mesh is refined. We see that the effectivities remain bounded between $$1$$ and $$20$$. This is a larger interval than in both Examples 1 and 2, although the largest effectivity occurs for $$1e-1$$. Fig. 3. View largeDownload slide Example 3. Interior layer problem under adaptive refinement. (a) Effectivities for varying $$\varepsilon$$. (b) Error convergence under uniform and adaptive refinement. Fig. 3. View largeDownload slide Example 3. Interior layer problem under adaptive refinement. (a) Effectivities for varying $$\varepsilon$$. (b) Error convergence under uniform and adaptive refinement. A comparison between uniform and adaptive refinement is shown in Fig. 3(b) for $$\varepsilon = 1e-2$$ and $$\varepsilon = 1e-4$$. In both cases, on later meshes, the adaptive refinement algorithm shows around 1 order of magnitude improvement in error for comparable degrees of freedom. 6.4 Comparison with other a posteriori bounds A comparison with other a posteriori estimates is in place. The pioneering works Sangalli (2001), Verfürth (2005), Sangalli (2008), Schötzau & Zhu (2009) introduce augmented norms for which a posteriori estimators are shown to admit so-called robust upper and lower bounds. In the context of dG methods, the augmented norm \[ | \kern -.25mm \|{w}| \kern -.25mm \|^2_{\rm aug} = \varepsilon| \kern -.25mm \|{w}| \kern -.25mm \|^2+\|cw\|^2 +\sum_{i=0}^N\left(ch+\frac{h}{\varepsilon}\right)[ \kern -.7mm [{w}] \kern -.7mm ]^2(x_i)+|bw|_{\ast}^2, \] where $$|q|_{\ast} := \sup_{v \in H^1_0(\alpha,\beta)\backslash\{0\}}\big(\int_{\alpha}^{\beta} qv' \,\mathrm{d} x)/(\varepsilon| \kern -.25mm \|{v}| \kern -.25mm \|^2+\|cv\|^2\big)^{1/2}$$ for $$q \in L^2(\alpha,\beta)$$ was used in Schötzau & Zhu (2009) to show that the estimator \begin{align}\label{sz} \eta^2_{\rm SZ} & = \sum_{i=1}^N\|c_{\rm osc}(\pi f+\varepsilon u_h^{\prime\prime}-bu^\prime_h-cu_h)\|^2_{T_i}+\sum_{i=1}^{N-1}\varepsilon^{-1/2}c_{\rm osc}[ \kern -.7mm [{u^\prime_h}] \kern -.7mm ]^2(x_i)\nonumber\\ &\quad{} +\sum_{i=0}^N\left(\varepsilon \sigma+ch+\frac{h}{\varepsilon}\right)[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i) \end{align} (6.1) is robust for $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\rm aug}^2$$, up to data oscillation; here we have taken $$c\ge 0$$ constant, for simplicity of the exposition. The inclusion of a dual-type $$|\cdot|_*$$-seminorm is an essential feature for the robustness of the estimator (6.1); corresponding observations are also true in the context of stabilized conforming methods (Sangalli, 2001, 2008; Verfürth, 2005). Indeed, going back to the setting of Example 1, the respective effectivities $${\eta}_{\rm SZ}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ with $${\eta}_{\rm SZ}$$ as in (6.1) and $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ the dG-energy norm defined in (2.5) (which, crucially, does not include the $$|\cdot|_*$$-seminorm) grow as $$\varepsilon$$ becomes smaller in the pre-asymptotic regime. In Fig. 4(a), we give these effectivities against the number of degrees of freedom $$N$$. Figures 4(b,c) show $$\eta_{\rm R}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ and $$\eta_{\rm J}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ against $$N$$, respectively, with \[ \eta^2_{\rm R}=\sum_{i=1}^Nc_{\rm osc}\|\pi f+\varepsilon u_h^{\prime\prime}-bu^\prime_h-cu_h\|^2_{T_i}\qquad\text{ and }\qquad \eta^2_{\rm J}=\sum_{i=0}^{N}\frac{h}{\varepsilon}[ \kern -.7mm [{u_h}] \kern -.7mm ]^2(x_i). \] Fig. 4. View largeDownload slide Example 1. Effectivity of the a posteriori estimator (6.1) under uniform refinement without the inclusion of the $$|\cdot|_*$$ in the error norm. (a) Effectivity of the $$\eta_{SZ}$$ estimator. (b) Contribution of the residual term $$\eta_{\rm R}$$. (c) Contribution of the jump term $$\eta_{\rm J}$$. Fig. 4. View largeDownload slide Example 1. Effectivity of the a posteriori estimator (6.1) under uniform refinement without the inclusion of the $$|\cdot|_*$$ in the error norm. (a) Effectivity of the $$\eta_{SZ}$$ estimator. (b) Contribution of the residual term $$\eta_{\rm R}$$. (c) Contribution of the jump term $$\eta_{\rm J}$$. Crucially, the estimator from Theorems 4.2 and 4.3 does not contain the terms $$\eta_{\rm R}$$ and $$\eta_{\rm J}$$ and, hence, robustness with respect to the standard dG-energy norm $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ is possible. To assess the potential impact of the above discussion in the context of adaptivity, we compare using $$\eta_{\rm SZ}$$ or $$\eta$$ to drive the adaptive algorithm for the problem of Example 2. A comparison of using $$\eta_{\rm SZ}$$ or $$\eta$$ to drive the adaptive algorithm is carried out for various values of $$\varepsilon$$ with the error measured in the dG-energy norm $$| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$; in Fig. 5, we show the results for $$\varepsilon=1e-5,1e-6,1e-7$$. We notice that $$\eta$$ initially leads to meshes that give a reduced error for a comparable number of degrees of freedom; indeed, the improvement is more marked for smaller $$\varepsilon$$. Fig. 5. View largeDownload slide Example 2. Performance of adaptive algorithm driven by $$\eta$$ and by $$\eta_{\rm SZ}$$. (a) $$\varepsilon = 1e-5$$. (b) $$\varepsilon = 1e-6$$. (c) $$\varepsilon = 1e-7$$. Fig. 5. View largeDownload slide Example 2. Performance of adaptive algorithm driven by $$\eta$$ and by $$\eta_{\rm SZ}$$. (a) $$\varepsilon = 1e-5$$. (b) $$\varepsilon = 1e-6$$. (c) $$\varepsilon = 1e-7$$. We also consider the estimator based on reconstructions of the diffusive and convective fluxes onto the Raviart–Thomas finite element space and on reconstructions of the potential into $$H^1_0(\it\Omega)$$, presented in Ern et al. (2010). Let the restriction of this estimator to the one-dimensional setting be denoted by $$\eta_{\rm ESV}$$; we omit full details of $$\eta_{\rm ESV}$$ for conciseness and we refer to Ern et al. (2010, Theorem 3.5). We investigate numerically whether $$\eta_{\rm ESV}$$ is also a robust estimator for the dG-energy norm in the present one-dimensional setting for the problem of Example 1. Figure 6 shows the effectivities $${\eta}_{\rm ESV}/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_{\gamma}$$ for $$\epsilon = 1\times 10^n$$, $$n=-4,\ldots,0$$, as the mesh is uniformly refined. We notice that, as $$\epsilon$$ is reduced, the maximum effectivity increases as in the case of $$\eta_{SZ}$$. Fig. 6. View largeDownload slide Effectivity of estimator under uniform refinement. (a) Effectivity of $$\eta_{\rm ESV}$$ under uniform refinement. (b) Effectivity of the element residual term $$\eta_{R,T}$$ (cf. Ern et al., 2010, Equation (35) for precise definition) of $$\eta_{\rm ESV}$$ under uniform refinement. Fig. 6. View largeDownload slide Effectivity of estimator under uniform refinement. (a) Effectivity of $$\eta_{\rm ESV}$$ under uniform refinement. (b) Effectivity of the element residual term $$\eta_{R,T}$$ (cf. Ern et al., 2010, Equation (35) for precise definition) of $$\eta_{\rm ESV}$$ under uniform refinement. 6.5 Example 4 To highlight the potential generality of the proposed estimator for higher-dimensional problems, we consider an error estimate analogous to that of Theorem 4.2 in the two-dimensional setting. As mentioned in the introduction, such a result can be obtained, provided we have at our disposal a technique that leads to error control of the limiting hyperbolic problem. Such results are available in certain cases (Georgoulis et al., 2014). To this end, we consider the problem \begin{eqnarray} -\varepsilon{\it\Delta} u+{\bf b}\cdot \nabla u + u &=& 1, \quad (x,y) \in {\it\Omega} = (0,1)^2, \nonumber\\ u &=& 0, \quad (x,y) \in \partial {\it\Omega}, \nonumber \end{eqnarray} where $${\bf b} = (\sqrt{2},1)^{\rm T}$$ and $$\varepsilon > 0$$; when $$\varepsilon=0$$, the boundary conditions are applied only on the inflow part of $$\partial{\it\Omega}$$. For $$\varepsilon > 0$$ boundary layers are formed along $$x=1$$ and $$y=1$$. We compare the error estimator under an adaptive refinement algorithm, for $$\varepsilon = 1e-2, 1e-4, 0$$, where a bulk marking strategy is applied; in all cases, an initial uniform grid comprising $$512$$ simplices is used. The effectivities $$\eta/| \kern -.25mm \|{u-u_h}| \kern -.25mm \|_1$$ are shown in Fig. 7, where we notice that in all cases they remain bounded between $$1$$ and $$4$$: for $$\varepsilon = 1e-2$$, the effectivities are slowly increasing but would appear to be converging to an upper bound below $$4$$, while for $$\varepsilon = 0$$ the effectivities seem to converge to $$1$$. The results are remarkably consistent with those from Example 2. In Fig. 8, adaptively refined grids are shown for both $$\varepsilon = 1e-2$$ and $$\varepsilon = 1e-4$$. We notice that for $$\varepsilon = 1e-2$$ refinement has been carried out in the vicinity of the boundary layers but also along a line parallel to the convection $$\bf b$$; for $$\varepsilon = 1e-4$$, where the boundary layers are narrower, refinement has been carried out only at the boundaries. Fig. 7. View largeDownload slide Example 4: (a) Error bound effectivities for varying $$\varepsilon$$ Fig. 7. View largeDownload slide Example 4: (a) Error bound effectivities for varying $$\varepsilon$$ Fig. 8. View largeDownload slide Example 4: Adaptive mesh for $$\varepsilon = 1e-2$$, iteration number 13 (left) and for $$\varepsilon = 1e-4$$, iteration number 20 (right). Fig. 8. View largeDownload slide Example 4: Adaptive mesh for $$\varepsilon = 1e-2$$, iteration number 13 (left) and for $$\varepsilon = 1e-4$$, iteration number 20 (right). 6.6 Example 5 In our final numerical experiment, we show that the estimator can drive an adaptive strategy to resolve solutions with both interior and boundary conditions. To this end, we consider the same test problem as in Schötzau & Zhu (2009, Example 3). In this case, the problem is to find $$u$$ such that \begin{alignat*}{2} -\varepsilon{\it\Delta} u+{\bf b}\cdot \nabla u &= 0, &(x,y) \in {\it \Omega} = (-1,1)^2, \\ u &= 0 &\text{on } x=-1 \text{ and } y = 1, \\ u &= \tanh\left(\frac{1-y}{\epsilon}\right) &\text{on $x=1$}, \\ u &= \frac{1}{2}\left(\tanh\left(\frac{x}{\epsilon}\right)+1\right) &\text{on $y=-1$}, \end{alignat*} where $${\bf b} = (-\sin\frac{\pi}{6},\cos\frac{\pi}{6})^{\rm T}$$. The boundary conditions cause $$u$$ to have an internal layer along the line $$y+\sqrt{3}x = -1$$ and boundary layers at the outflow boundaries. Since there is no exact solution available for this problem, reference solutions for varying $$\varepsilon$$ have been obtained using a different a posteriori estimator (namely, the $$\eta_{\rm SZ}$$ estimator) within an adaptive framework and setting $$p=4$$ to ensure that the reference solution is of sufficient accuracy. In this case, we apply a fixed fraction strategy with $$25\%$$ of the elements chosen for refinement and $$10\%$$ chosen for derefinement at each iteration. For both $$\varepsilon = 10^{-2}$$ and $$\varepsilon = 10^{-4}$$, we begin with a uniform mesh comprising $$512$$ right-angled triangles, which are not aligned with the direction of flow, and set $$p=1$$. Figure 9 reports the error convergence history for both these cases, while Fig. 10 shows the resultant meshes after $$14$$ iterations of the adaptive algorithm for $$\varepsilon = 10^{-2}$$ and $$\varepsilon = 10^{-4}$$, respectively. Fig. 9. View largeDownload slide Example 5: Error convergence for adaptive strategy for $$\varepsilon=1e-2$$ and $$\varepsilon=1e-4$$. Fig. 9. View largeDownload slide Example 5: Error convergence for adaptive strategy for $$\varepsilon=1e-2$$ and $$\varepsilon=1e-4$$. Fig. 10. View largeDownload slide Example 5: Adaptive mesh for $$\varepsilon = 1e-2$$ (left) and for $$\varepsilon = 1e-4$$ (right) at iteration number $$14$$. Fig. 10. View largeDownload slide Example 5: Adaptive mesh for $$\varepsilon = 1e-2$$ (left) and for $$\varepsilon = 1e-4$$ (right) at iteration number $$14$$. The error indicator leads to refinement in the regions of both the boundary and interior layers for both values of $$\varepsilon$$. For $$\varepsilon=10^{-4}$$, a pre-asymptotic regime is observed, as expected; when the mesh is fine enough to resolve the boundary and interior layers, the error begins to converge at the same rate as for $$\varepsilon = 10^{-2}$$. References Ainsworth M. , Allendes A. , Barrenechea G. R. & Rankin R. ( 2013 ) Fully computable a posteriori error bounds for stabilized FEM approximations of convection-reaction-diffusion problems in three dimensions. Internat. J. Numer. Methods Fluids 73 , 765 – 790 . Arnold D. N. ( 1982 ) An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal. , 19 , 742 – 760 . Google Scholar CrossRef Search ADS Cohen A. , Dahmen W. & Welper G. ( 2012 ) Adaptivity and variational stabilization for convection-diffusion equations. ESAIM Math. Model. Numer. Anal. , 46 , 1247 – 1273 . Google Scholar CrossRef Search ADS Creusé E. & Nicaise S. ( 2013 ) A posteriori error estimator based on gradient recovery by averaging for convection-diffusion-reaction problems approximated by discontinuous Galerkin methods. IMA J. Numer. Anal. , 33 , 212 – 241 . Google Scholar CrossRef Search ADS Dahmen W. , Huang C. , Schwab C. & Welper G. ( 2012 ) Adaptive Petrov-Galerkin methods for first order transport equations. SIAM J. Numer. Anal. , 50 , 2420 – 2445 . Google Scholar CrossRef Search ADS Ern A. , Stephansen A. F. & Vohralík M. ( 2010 ) Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems. J. Comput. Appl. Math. , 234 , 114 – 130 . Google Scholar CrossRef Search ADS Georgoulis E. H. , Hall E. & Makridakis C. ( 2014 ) Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems . Barrett Lectures on Discontinuous Galerkin Methods. IMA Vol. Math. Appl. , vol. 157 , Cham : Springer , 195 – 207 . Georgoulis E. H. & Lasis A. ( 2006 ) A note on the design of $$hp$$-version interior penalty discontinuous Galerkin finite element methods for degenerate problems. IMA J. Numer. Anal. , 26 , 381 – 390 . Google Scholar CrossRef Search ADS Houston P. , Schwab C. & Süli E. ( 2002 ) Discontinuous $$hp$$-finite element methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. , 39 , 2133 – 2163 (electronic). Google Scholar CrossRef Search ADS Johnson C. & Pitkäranta J. ( 1986 ) An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math. Comput. , 46 , 1 – 26 . Google Scholar CrossRef Search ADS Makridakis C. & Nochetto R. H. ( 2006 ) A posteriori error analysis for higher order dissipative methods for evolution problems. Numer. Math. , 104 , 489 – 514 . Google Scholar CrossRef Search ADS Roos H.-G. , Stynes M. & Tobiska L. ( 2008 ) Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems , vol. 24 , 2nd edn . Springer Series in Computational Mathematics , Berlin : Springer . Sangalli G. ( 2001 ) A robust a posteriori estimator for the residual-free bubbles method applied to advection-diffusion problems. Numer. Math. , 89 , 379 – 399 . Google Scholar CrossRef Search ADS Sangalli G. ( 2008 ) Robust a-posteriori estimator for advection-diffusion-reaction problems. Math. Comp. , 77 , 41 – 70 (electronic). Google Scholar CrossRef Search ADS Schötzau D. & Wihler T. P. ( 2010 ) A posteriori error estimation for $$hp$$-version time-stepping methods for parabolic partial differential equations. Numer. Math. , 115 , 475 – 509 . Google Scholar CrossRef Search ADS Schötzau D. & Zhu L. ( 2009 ) A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations. Appl. Numer. Math. , 59 , 2236 – 2255 . Google Scholar CrossRef Search ADS Verfürth R. ( 2005 ) Robust a posteriori error estimates for stationary convection-diffusion equations. SIAM J. Numer. Anal. , 43 , 1766 – 1782 (electronic). Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Dec 22, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off