Eliminating Gibbs Phenomena: A Non-linear Petrov-Galerkin Method for the Convection-Diffusion-Reaction Equation

TL;DR

This paper introduces a non-linear Petrov-Galerkin method based on residual minimization in dual norms within Banach spaces to effectively eliminate Gibbs phenomena in convection-diffusion-reaction equations, especially in convection-dominated regimes.

Contribution

It extends the Petrov-Galerkin framework to reflexive Banach spaces with residual minimization in dual norms, reducing oscillations in numerical solutions of convection-diffusion-reaction equations.

Findings

Oscillations diminish as q approaches 1 in the Lq setting.

The method effectively reduces Gibbs phenomena in numerical experiments.

Generalizes discontinuous Petrov-Galerkin methods to Banach spaces.

Abstract

In this article we consider the numerical approximation of the convection-diffusion-reaction equation. One of the main challenges of designing a numerical method for this problem is that boundary layers occurring in the convection-dominated case can lead to non-physical oscillations in the numerical approximation, often referred to as Gibbs phenomena. The idea of this article is to consider the approximation problem as a residual minimization in dual norms in Lq-type Sobolev spaces, with 1 < q < $\infty$. We then apply a non-standard, non-linear PetrovGalerkin discretization, that is applicable to reflexive Banach spaces such that the space itself and its dual are strictly convex. Similar to discontinuous Petrov-Galerkin methods, this method is based on minimizing the residual in a dual norm. Replacing the intractable dual norm by a suitable discrete dual norm gives rise to a non-linear inexact mixed method. This generalizes the Petrov-Galerkin framework developed in the context of discontinuous Petrov-Galerkin methods to more general Banach spaces. For the convection-diffusion-reaction equation, this yields a generalization of a similar approach from the L2-setting to the Lq-setting. A key advantage of considering a more general Banach space setting is that, in certain cases, the oscillations in the numerical approximation vanish as q tends to 1, as we will demonstrate using a few simple numerical examples.

Equations

$$⟨\phi ,v{⟩}_{V^{\prime },V}:=\phi (v).$$

$$\Vert \phi {\Vert }_{V^{\prime }}:=\underset{v\in V,v\ne 0}{\sup }\frac{⟨\phi ,v{⟩}_{V^{\prime },V}}{\Vert v{\Vert }_V}.$$

$$Bu=f\text{ in }{V}^{\prime },$$

$$u_n=\arg \underset{w_n\in U_n}{\min }\Vert f-B{w}_n{\Vert }_{V^{\prime }}.$$

$$\Vert u-u_n{\Vert }_U\le \frac{M}{\gamma }\underset{w_n\in U_n}{\inf }\Vert u-w_n{\Vert }_U,$$

$$\begin{array}{rlrl}⟨J_V(r),v{⟩}_{V^{\prime },V}+⟨B{u}_n,v{⟩}_{V^{\prime },V}& =⟨f,v{⟩}_{V^{\prime },V},& & \forall v\in V,\\ ⟨B{w}_n,r{⟩}_{V^{\prime },V}& =0,& & \forall w_n\in U_n,\end{array}$$

$$\begin{array}{rl}-\epsilon \Delta u+\mathbf{b}\cdot \nabla u+cu& =f\text{ in }\Omega ,\\ u& =0\text{ on }\Gamma =\partial \Omega ,\end{array}$$

$${\mathcal{B}}_{\epsilon }(u,v)$$

$${\mathcal{J}}_V^{\phi }(v):=\left\{v^{\prime }\in V^{\prime }:⟨v^{\prime },v{⟩}_{V^{\prime },V}=\Vert v{\Vert }_V\Vert v^{\prime }{\Vert }_{V^{\prime }},\Vert v^{\prime }{\Vert }_{V^{\prime }}=\phi (\Vert v{\Vert }_V)\right\}$$

$${\mathcal{J}}_V^{\phi }(v)=\partial F_V^{\phi }(v),$$

$$\Vert v{\Vert }_{L^q\left(\Omega \right)}={\left({\int }_{\Omega }|v{|}^q\mathrm{d}\mathbf{x}\right)}^{1/q}.$$

$$\begin{array}{rl}⟨J_q(v),w{⟩}_{L^{q^{\prime }}\left(\Omega \right),L^q\left(\Omega \right)}& ={\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{q}\Vert v+tw{\Vert }_{L^q\left(\Omega \right)}^q\right)|}_{t=0}\\ & ={\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{q}{\int }_{\Omega }|v+tw{|}^q\mathrm{d}\mathbf{x}\right)|}_{t=0}\\ & ={\int }_{\Omega }|v{|}^{q-1}\mathrm{sgn}(v)w\mathrm{d}\mathbf{x}.\end{array}$$

$$\Vert J_q(v){\Vert }_{L^{q^{\prime }}\left(\Omega \right)}=\Vert v{\Vert }_{L^q\left(\Omega \right)}^{q-1},⟨J_q(v),v{⟩}_{L^{q^{\prime }}\left(\Omega \right),L^q\left(\Omega \right)}=\Vert v{\Vert }_{L^q\left(\Omega \right)}^q.$$

$$⟨{\tilde{J}}_q(v),w{⟩}_{L^{q^{\prime }}\left(\Omega \right),L^q\left(\Omega \right)}$$

$$⟨{\tilde{J}}_q^{\prime }(v)(z),w{⟩}_{L^{q^{\prime }}\left(\Omega \right),L^q\left(\Omega \right)}$$

$$+(2-q)\Vert v{\Vert }_{L^q\left(\Omega \right)}^{2-2q}{\int }_{\Omega }|v{|}^{q-2}vz\mathrm{d}\mathbf{x}{\int }_{\Omega }|v{|}^{q-2}vw\mathrm{d}\mathbf{x}.$$

$$⟨J_{W_0^{1,q}\left(\Omega \right)}(v),w{⟩}_{L^{q^{\prime }}\left(\Omega \right),L^q\left(\Omega \right)}={\int }_{\Omega }\sum _{i=1}^d|{\partial }_{i}v{|}^{q-1}\mathrm{sgn}({\partial }_{i}v){\partial }_{i}w\mathrm{d}\mathbf{x}.$$

$$\Vert x{\Vert }_{l^q}$$

$$\Vert x{\Vert }_{l^r}$$

$${\left({\int }_{\Omega }\Vert \nabla v{\Vert }_{l^2}^q\mathrm{d}\mathbf{x}\right)}^{1/q}$$

$${\int }_{\Omega }\Vert \nabla v{\Vert }_{l^2}^{q-2}\nabla v\cdot \nabla w\mathrm{d}\mathbf{x},$$

$$\underset{w\in U}{\inf }\underset{v\in V}{\sup }\frac{b(w,v)}{\Vert w{\Vert }_U\Vert v{\Vert }_V}=\gamma >0,$$

$$\left\{v\in V:b(w,v)=0,\forall w\in U\right\}=\{0\}.$$

$$b(u,v)=\ell (v)\forall v\in V.$$

$$\begin{array}{rl}{u}_n& =\arg \underset{w_n\in U_n}{\min }{H}_V(w_n),\\ H_V(w_n)& :=\psi (\Vert B{w}_n-\ell {\Vert }_{V^{\prime }}),\end{array}$$

$$\begin{array}{rl}⟨J_{V^{\prime }}^{\phi }(\ell -B{u}_n),& B{w}_n{⟩}_{V^{\prime \prime },V^{\prime }}\\ & ={⟨B{w}_n,{\left(J_V^{{\phi }^{-1}}\right)}^{-1}(\ell -B{u}_n)⟩}_{V^{\prime },V}=0\forall w_n\in U_n.\end{array}$$

$$⟨J_V^{{\phi }^{-1}}(r),v{⟩}_{V^{\prime },V}+⟨B{u}_n,v{⟩}_{V^{\prime },V}$$

$$⟨B{w}_n,r{⟩}_{V^{\prime },V}$$

$$\Vert u-u_n{\Vert }_U\le \frac{1}{\gamma }\underset{v\in V}{\sup }\frac{b(u-u_n,v)}{\Vert v{\Vert }_V}=\frac{1}{\gamma }\Vert \ell -B{u}_n{\Vert }_{V^{\prime }}\le \frac{M}{\gamma }\Vert u-u_n{\Vert }_U,$$

$$\Vert v{\Vert }_{\mathrm{opt}}:=\underset{u\in U}{\sup }\frac{b(u,v)}{\Vert u{\Vert }_U}=\Vert B^{\prime }v{\Vert }_{U^{\prime }}.$$

$$\underset{v\in V}{\sup }\frac{b(w,v)}{\Vert v{\Vert }_{\mathrm{opt}}}=\underset{v\in V}{\sup }\frac{(B^{\prime }v)(w)}{\Vert B^{\prime }v{\Vert }_{U^{\prime }}}=\underset{g\in U^{\prime }}{\sup }\frac{g(w)}{\Vert g{\Vert }_{U^{\prime }}}=\Vert w{\Vert }_U.$$

$$\Vert u{\Vert }_E:=\underset{v\in V}{\sup }\frac{b(u,v)}{\Vert v{\Vert }_V}.$$

Full text

Keywords: convection-diffusion; Petrov-Galerkin; Gibbs phenomenon; finite element methods; Banach spaces

AMS Subject Classification: 65N30; 35J20

1 Introduction

The term Gibbs phenomenon originally refers to the effect that the partial sums of the Fourier series approximating a function with a jump discontinuity exhibit over- and undershoots near the discontinuity. The phenomenon is named after Willard Gibbs who described it in 1899 [1], though it had already been discovered earlier by Henry Wilbraham in 1848 [2]. It also occurs in the best approximation by spline functions in the $L^{2}$-metric [3] and is one of the main challenges in the numerical approximation of partial differential equations (PDEs) whose solutions contain sharp features such as shocks or thin layers. In [4] it is shown in one dimension that the best approximation of jump discontinuities by polygonal lines on a uniform grid in one dimension does not lead to Gibbs phenomenon in $L^{1}$. More precisely, it is shown that the overshoot of the best approximation in $L^{q}$, $1<q<\infty$, is an increasing function of $q$ that tends to [math] as $q$ tends to $1$. We show in [5] that this result is still true for the best approximation in $L^{q}$ by piecewise linear functions on certain meshes in two dimensions and certain non-uniform meshes in one dimension. However, there exist meshes both in one and two dimensions such that the (maximal) overshoot tends to some $\alpha>0$ as $q\rightarrow 1$. Nonetheless, it is suggested in [5] that if the location of the discontinuity is known, a mesh can be constructed in such a way that the overshoot vanishes. In [6], it is shown that Gibbs phenomena occur in the best approximation of a function containing a jump discontinuity by a trigonometric polynomial in the $L^{1}$-metric. This confirms that the choice of the finite dimensional approximation space crucially determines whether Gibbs phenomena can be eliminated by considering best approximation in $L^{1}$.

The motivation of this article is to exploit the Gibbs-reducing property of $L^{q}$-type spaces in the numerical approximation of PDEs using finite element methods. A similar idea has already been pursued for numerical methods in $L^{q}$ by Guermond [7]. In [7], Guermond points out that there are only very few attempts to approximate PDEs directly in $L^{1}$ despite the fact that first-order PDEs and their non-linear generalizations have been extensively studied in $L^{1}$. The existing numerical methods to achieve this include the ones outlined in the articles by Lavery [8, 9, 10], the reweighted least-squares method of Jiang [11, 12] and the methods outlined in the series of articles by Guermond et al. [7, 13, 14, 15, 16].

More recently, a novel approach to designing finite element methods in a very general Banach space setting has been introduced in [17]. This approach is rooted in the so-called Discontinuous Petrov-Galerkin (DPG) methods [18] and extends the concept of optimal test norms and functions from Hilbert spaces to more general Banach spaces yielding a scheme that can be interpreted as a non-linear Petrov-Galerkin method. At least in an abstract sense, this approach outlines how to design a numerical method that leads to a quasi-best approximation of the solution in a space of choice, provided the continuous problem is well-posed in a suitable sense. In this article we apply this abstract approach to the convection-diffusion-reaction equation in a $W^{1,q}_{0}\left(\Omega\right)$-$W^{1,q^{\prime}}_{0}\left(\Omega\right)$ setting where $1/q+1/q^{\prime}=1$ and study the effect this has on the numerical approximation of boundary layers. We will see that as $q\rightarrow 1$, the Gibbs phenomenon can be eliminated entirely on some meshes. We will also consider certain choices of meshes for which this is not the case and consider how this can be fixed.

1.1 Notation

Throughout this article, we denote by $L^{q}\left(\Omega\right)$, $1\leq q<\infty$, the Lebesgue space of $q$-integrable functions on a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$, $d\in\{1,2,3,\ldots\}$; $L^{\infty}(\Omega)$ is the Lebesgue space of functions on $\Omega$ with finite essential supremum; and $W^{1,q}\left(\Omega\right)$, $1\leq q\leq\infty$, is the Sobolev space of functions that are in $L^{q}\left(\Omega\right)$ such that their gradient is in $L^{q}\left(\Omega\right)^{d}$. Furthermore, $W^{1,q}_{0}\left(\Omega\right)\subset W^{1,q}\left(\Omega\right)$ is the subspace of all functions with zero trace on the boundary $\partial\Omega$. The corresponding norms are denoted by $\|\cdot\|_{L^{q}\left(\Omega\right)}$ and $\|\cdot\|_{W^{1,q}\left(\Omega\right)}$, respectively, and the Sobolev-semi norm on $W^{1,q}\left(\Omega\right)$ is given by $|\cdot|_{W^{1,q}\left(\Omega\right)}$. For $q=2$, we furthermore use the usual notation $H^{1}(\Omega):=W^{1,2}\left(\Omega\right)$ and $H^{1}_{0}\left(\Omega\right):=W^{1,2}_{0}\left(\Omega\right)$. For $1\leq q\leq\infty$ we write $q^{\prime}$ to denote the dual exponent such that $1/q+1/q^{\prime}=1$. For any Banach space $V$, its norm is denoted by $\|\cdot\|_{V}$ and its dual space by $V^{\prime}$; the dual space of $W^{1,q}_{0}\left(\Omega\right)$ is given by $W^{-1,q^{\prime}}\left(\Omega\right)$ and $H^{-1}(\Omega):=W^{-1,2}\left(\Omega\right)$. For $v\in V$ and $\varphi\in V^{\prime}$, we have the duality pairing

[TABLE]

For any $\varphi\in V^{\prime}$, its norm in the dual space $V^{\prime}$ is given by

[TABLE]

If $V$ is a Hilbert space, we denote the Riesz map on $V$ by $R_{V}$.

1.2 Problem Statement

The approach to generalizing the DPG framework to general Banach spaces introduced in [17] is based on residual minimization. To this end, the authors consider the abstract problem: find $u\in U$ such that

[TABLE]

where $U$ and $V$ are Banach spaces, $B:U\rightarrow V^{\prime}$ is a continuous, bounded-below linear operator, and the right hand side $f$ is in the dual space $V^{\prime}$. The associated residual minimization problem for a given finite dimension subspace $U_{n}\subset U$, e.g., a finite element space, of dimension $n$ is defined as follows:

[TABLE]

It can be shown that the solution $u_{n}$ is a quasi-best approximation of the analytical solution $u$, in the sense that

[TABLE]

where $M$ is the continuity constant of $B$ and $\gamma$ the bounded-below constant of $B$.

In [17] it is shown that the minimization problem can be reformulated yielding the following non-linear mixed system: find $(r,u_{n})\in V\times U_{n}$ such that

[TABLE]

where $J_{V}$ is a so-called duality mapping. Duality mappings are a generalization of the Riesz map $R_{V}$ to general Banach spaces with the main difference that they are non-linear mappings unless $V$ is a Hilbert space. To turn the above mixed system into a practical method, $V$ is additionally replaced by a finite dimensional subspace $V_{m}$ of dimension $m$. As a result the minimization problem is no longer solved exactly. However, if the spaces $U_{n}$ and $V_{m}$ are chosen in a suitable way, one can obtain a well-posed fully discrete mixed system that retains the quasi-best approximation property of $u_{n}$ with a modified constant.

The aim of this article is to apply this abstract framework to the convection-diffusion-reaction equation

[TABLE]

where $\varepsilon$, $b:\Omega\rightarrow\mathbb{R}^{d}$, and $c:\Omega\rightarrow\mathbb{R}$ are the (positive) diffusion parameter, convection field and reaction coefficient, respectively, and $f:\Omega\rightarrow\mathbb{R}$ is a given source term. Therefore, we will consider the linear operator associated with the following bilinear form:

[TABLE]

The resulting method is a generalization to the $W^{1,q}_{0}\left(\Omega\right)$-$W^{1,q^{\prime}}_{0}\left(\Omega\right)$-setting of the method introduced for the $H^{1}_{0}\left(\Omega\right)$-setting in [19], where the DPG framework is applied to the convection-diffusion-reaction equation without introducing broken test spaces. Related Petrov-Galerkin formulations are studied in [20, 21, 22, 23, 24].

1.3 Summary of Results

In this article we study how the mixed method (1.1) and its fully discrete counterpart can be applied to the convection-diffusion-reaction equation (1.2). In particular, we study the choice of the space $V$, the corresponding norm $\|\cdot\|_{V}$, and the discrete space $V_{m}$. The norm $\|\cdot\|_{V}$ on the space $V$ determines the operator $J_{V}$ in (1.1) and thus crucially determines the resulting method. As in the special case $q=2$ investigated in [19], we introduce weakly imposed boundary conditions on the inflow boundary in the space $V_{m}$ to address robustness issues. We demonstrate in one dimension, how this choice of boundary conditions affects the approximation and compare this approach with the case when a weighted norm on $V$ is employed.

The main focus of our numerical investigation is eliminating Gibbs phenomena in the finite element approximation. Indeed, we will see that the numerical approximation generated by our method qualitatively behaves like the $L^{q}\left(\Omega\right)$-best approximation of the analytical solution. Thus, the Gibbs phenomenon can be eliminated by taking the limit $q\rightarrow 1$ provided that the $L^{1}\left(\Omega\right)$-best approximation does not exhibit Gibbs phenomena. The results in [5] show that this depends on the mesh that is chosen. In one dimension only certain non-uniform grids have the property that the $L^{1}\left(\Omega\right)$-best approximation contains over- and undershoots, whereas in higher dimensions this can even occur on structured, uniform meshes. We will demonstrate for one two-dimensional example that it is possible to use the insights from [5] to modify the mesh in order to eliminate Gibbs phenomena as $q\rightarrow 1$.

1.4 Outline of the Paper

This article is structured as follows: in Section 2 we introduce duality mappings, which are essential to the abstract framework in [17], cf., (1.1). In Section 3 we give a brief overview of the abstract framework in [17]. Then, in Section 4 we then apply this framework to the convection-diffusion-reaction equation and illustrate how this yields a practical method that can be implemented, e.g., in FEniCS [25, 26]. In Section 5 we investigate the practical performance of the proposed method for a range of different test problems in one and two dimensions. Finally, in 6, we summarize the work presented in this article and highlight the potential and challenges of the proposed numerical method.

2 Duality Mappings

One of the key ingredients used to extend concepts from Hilbert spaces to more general Banach spaces both in [27] and [28] is replacing the Riesz map with so-called duality mappings. In general, duality mappings are non-linear maps that share certain properties with the Riesz map. If the underlying space is a Hilbert space, any duality mapping is linear and the Riesz map is one possible choice for the duality mapping. In this section we will give a precise definition of duality mappings, summarize a range of very useful properties and list some relevant examples.

Definition 2.1**.**

A weight function is a continuous and strictly increasing function $\varphi:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ such that $\varphi(0)=0$ and $\lim_{t\rightarrow\infty}\varphi(t)=+\infty$. 2. 2.

Let $V$ be a Banach space and $\varphi$ a weight function. Denote by $\mathcal{P}({V}^{\prime})$ the power set of ${V}^{\prime}$. Then the multivalued map $\mathcal{J}_{V}^{\varphi}:V\rightarrow\mathcal{P}({V}^{\prime})$, defined by

[TABLE]

is called a duality mapping of weight $\varphi$.

Due to a corollary of the Hahn-Banach theorem (cf., e.g., [29, Corollary 1.3]) the set $\mathcal{J}^{\varphi}_{V}(v)$ is non-empty. If we choose $\varphi(t)=t$ as a weight function we obtain the so-called normalized duality map. In many cases this choice is the simplest and most convenient one. In fact, some books only introduce the duality mapping for this choice of $\varphi$, cf., e.g., [29, 30, 31]. As we will see later, however, $\varphi(t)=t^{q-1}$ is in some cases a more suitable choice for $L^{q}$-type spaces. If ${q}=2$, this choice coincides with the normalized duality map. Furthermore, note that we have $\mathcal{J}_{V}(v)=\{R_{V}(v)\}$ if $\mathcal{J}_{V}$ is the normalized duality map and $V$ is a Hilbert space due to the Riesz Representation theorem.

2.1 Properties of Duality Mappings

In the following proposition we summarize a few properties of the duality map for Banach spaces with particular structures.

Proposition 2.2**.**

Let $V$ be a Banach space and denote by $\mathcal{J}_{V}^{\varphi}:V\rightarrow\mathcal{P}({V}^{\prime})$ the duality map of weight $\varphi$ on $V$. Then the following statements are true:

${V}^{\prime}$* is strictly convex111A Banach space $V$ is strictly convex if for all $v_{1},v_{2}\in V$ such that $v_{1}\neq v_{2}$ and $\|v_{1}\|_{V}=\|v_{2}\|_{V}=1$ it holds that $\|\vartheta v_{1}+(1-\vartheta)v_{2}\|_{V}<1\quad\forall\,\vartheta\in(0,1).$ if and only if $\mathcal{J}_{V}^{\varphi}$ is single valued, cf., [30, Prop. 12.3]. In this case we define the duality map $J_{V}^{\varphi}:V\rightarrow{V}^{\prime}$ such that $\mathcal{J}_{V}^{\varphi}(v)=\{J_{V}^{\varphi}(v)\}$ for all $v\in V$.* 2. 2.

If $V$ is strictly convex, then $\mathcal{J}_{V}^{\varphi}(v)\cap\mathcal{J}_{V}^{\varphi}(w)=\emptyset$ for all $w\neq v$. In particular, $\mathcal{J}_{V}^{\varphi}$ is injective. 3. 3.

$V$* is reflexive if and only if $\mathcal{J}_{V}^{\varphi}$ is surjective in the sense that for every ${v}^{\prime}\in{V}^{\prime}$ there is a $v\in V$ such that ${v}^{\prime}\in\mathcal{J}_{V}^{\varphi}(v)$, cf., [32, Theorem 3.4, Chapter II].* 4. 4.

If $V$ is a reflexive Banach space and $\mathcal{J}^{\varphi}_{V}$ is a duality mapping of weight $\varphi$, then $(\mathcal{J}^{\varphi}_{V})^{-1}$ is a duality mapping on $V^{\prime}$ of weight $\varphi^{-1}$, cf., **[32, Cor. 3.5, Ch. II]**.

The main implication of the above proposition is that in a strictly convex and reflexive space $V$, the duality mapping is bijective and its inverse can be identified with a duality mapping on the dual space $V^{\prime}$ with the inverse weight $\varphi^{-1}$.

The following theorem is a special case of Theorem 4.4 in [32, Chapter I] and states that the duality map on $V$ can be characterized using the subdifferential of the norm on $V$. This is a key property of the duality map that will allow us to derive the duality map for some specific Banach spaces in the special case that the subdifferential is essentially the Gâteaux or Fréchet derivative of the norm.

Theorem 2.3** (Asplund, cf., [32, Ch. I, Theorem 4.4]).**

Let $V$ be a Banach space and define $F_{V}^{\varphi}:V\rightarrow\mathbb{R}$ by $F_{V}^{\varphi}(\cdot):=\psi(\|\cdot\|_{V})$, where $\psi(s):=\int_{0}^{s}\varphi(t)\mathop{}\!\mathrm{d}t$ and $\varphi$ is a weight function. Then for any $v\in V$, we have

[TABLE]

where $\partial F_{V}^{\varphi}(v)$ denotes the subdifferential of $F_{V}^{\varphi}$ at $v$.

The result that finally allows us to compute duality maps is the following.

Proposition 2.4** (cf., [31, Proposition 47.19]).**

Let $V$ be a Banach space and define $F_{V}^{\varphi}:V\rightarrow\mathbb{R}$ by $F_{V}^{\varphi}(\cdot):=\psi(\|\cdot\|_{V})$, where $\psi(s):=\int_{0}^{s}\varphi(t)\mathop{}\!\mathrm{d}t$ and $\varphi$ is a weight function.

If ${V}^{\prime}$ is strictly convex, then $\nabla F_{V}^{\varphi}(v)$ exists as a Gâteaux derivative and $\nabla F_{V}^{\varphi}(v)=J_{V}^{\varphi}(v)$ for all $v\in V$. 2. 2.

If ${V}^{\prime}$ is uniformly convex, then $\nabla F_{V}^{\varphi}(v)$ exists as a Fréchet derivative and $\nabla F_{V}^{\varphi}(v)=J_{V}^{\varphi}(v)$ for all $v\in V$.

Note that uniform convexity of ${V}^{\prime}$ implies strict convexity of ${V}^{\prime}$ (cf., [33]) and due to the Milman-Pettis Theorem, ${V}^{\prime}$ and hence also $V$ are reflexive in this case. Furthermore, if $F_{V}^{\varphi}$ is Gâteaux differentiable, Theorem 2.3 guarantees that the duality map is single valued which implies strict convexity of the dual space ${V}^{\prime}$ due to Proposition 2.2.

2.2 Some Examples of Duality Mappings on Sobolev Spaces

First, consider $V=L^{{q}}\left(\Omega\right)$ with the norm

[TABLE]

Let us denote the duality mapping on $L^{{q}}\left(\Omega\right)$ with weight $\varphi(t)=t^{{q}-1}$ by $J_{q}$. In this case $\psi(s)=\int_{0}^{s}\varphi(t)\mathop{}\!\mathrm{d}t=\frac{1}{{q}}s^{q}$ and thus we can compute

[TABLE]

Note that $\varphi^{-1}(t)=t^{{q^{\prime}}-1}$ for $q$ such that $1=1/q+1/q^{\prime}$ and therefore $J_{q^{\prime}}^{-1}=J_{q}$. Moreover, we have for all $v\in L^{q}\left(\Omega\right)$

[TABLE]

Similarly, we can compute the normalized duality map $\tilde{J}_{q}$ on $L^{{q}}\left(\Omega\right)$, i.e., the duality map with weight $\varphi(t)=t$, and obtain

[TABLE]

Comparing the expressions for $J_{q}$ and $\tilde{J}_{q}$, we can see that it may be useful, in particular for the implementation, to use $J_{q}$ instead of the normalized duality mapping in order to avoid the additional scaling $\|v\|_{L^{{q}}\left(\Omega\right)}^{2-{q}}$. Indeed, if $\tilde{J}_{q}$ is a non-linear operator in a variational problem that is approximated using a finite element method, we would need to evaluate the derivative of $\tilde{J}_{q}$ to determine the Jacobi matrix used in Newton’s method. The derivative is given by

[TABLE]

If $w$ and $z$ are finite element functions with local support, the first term is only non-zero if both $w$ and $z$ are non-zero and thus this term would yield a sparse matrix for typical finite element spaces. The second term on the other hand is always non-zero if $v$ is non-zero on the whole domain and thus may lead to a dense Jacobi matrix. If we consider $J_{q}$ instead, the derivative consists of only the first term without the scaling $\|v\|_{L^{q}\left(\Omega\right)}^{2-q}$ and therefore we obtain a sparse Jacobi matrix.

As a second example consider the space $W^{1,{q}}_{0}\left(\Omega\right)$ with the (semi-)norm $|v|_{W^{1,q}\left(\Omega\right)}$. In the same way as before we can compute the duality map of weight $\varphi(t)=t^{{q}-1}$ and obtain

[TABLE]

Remark 2.5* (q-Laplacian).*

Let $1\leq q<r<\infty$ and $x\in\mathbb{R}^{d}$. In this case we have the following norm equivalence on the finite dimensional space $\mathbb{R}^{d}$:

[TABLE]

As a result

[TABLE]

defines a norm on $W^{1,q}_{0}\left(\Omega\right)$ that is equivalent to $|v|_{W^{1,q}\left(\Omega\right)}$ and the corresponding duality mapping for $q>2$ is given by

[TABLE]

which is also the weak form of the $q$-Laplacian.

3 Minimum Residual Methods in Banach Spaces

As mentioned before, the method we are considering is a generalization of the scheme discussed in [19]. Conceptually, this method is closely related to DPG methods with the main difference that only continuous finite element spaces are considered in [19]. The key observation that makes it possible to extend the methodology to general Banach spaces is the equivalence with a residual minimization problem, cf., [19] and [34] for the DPG method. In this section, we will see that, in the context of Banach spaces, the Riesz map can be replaced by a duality mapping, but due to the non-linearity of the duality mapping we lose the concept of an optimal test space which is used in the context of DPG methods, cf., [35]. The extension to Banach spaces was first introduced in [27] and we will repeat the main concepts in this section.

Let $U$ and $V$ be two Banach spaces and $b:U\times V\rightarrow\mathbb{R}$ a continuous bilinear form that satisfies the following inf-sup conditions:

[TABLE]

For a given right-hand side $\ell\in V^{\prime}$, we consider the problem: find $u\in U$ such that

[TABLE]

Introducing a finite dimensional subspace $U_{n}\subset U$, we can formulate the following residual minimization problem: find $u_{n}\in U_{n}$ such that

[TABLE]

where $\psi(0)=0$, $\varphi:=\psi^{\prime}$ is a weight function as defined in Section 2 and $B:U\rightarrow V^{\prime}$ denotes the linear operator associated with the bilinear form $b$. A typical choice for $\psi$ would be $\psi(t)=t^{q}/q$, $1<q<\infty$. If ${V}^{\prime}$ is strictly convex, then the duality mapping $J_{{V}^{\prime}}^{\varphi}$ is single valued and $J_{{V}^{\prime}}^{\varphi}=\nabla(v^{\prime}\mapsto\psi(\|v^{\prime}\|_{{V}^{\prime}}))$.

Note that this best approximation problem can be formulated in any Banach space. However, we need strict convexity for uniqueness of minimizers, cf., e.g., [36, 27]. The closed range theorem is still applicable in general Banach spaces and provides conditions for well-posedness of linear problems. Its reformulation in terms of inf-sup conditions only requires reflexivity of the test space to describe the kernel of $B^{\prime}$ by (3.1b). Here, $B^{\prime}:V\rightarrow U^{\prime}$ denotes the adjoint operator to $B$.

3.1 Saddle Point Formulation

If ${V}^{\prime}$ is strictly convex, it can be shown that the residual minimization problem is equivalent to $\nabla H_{V}(u_{n})=J_{V^{\prime}}^{\varphi}(\ell-Bu_{n})=0$. If $V$ is additionally reflexive, as well as being strictly convex, we obtain the following non-linear Petrov Galerkin formulation: find $u_{n}\in U_{n}$ such that

[TABLE]

Here, we used $J^{\varphi}_{V^{\prime}}=\left(J_{{V}}^{\varphi^{-1}}\right)^{-1}$ by means of canonical identification, cf., Proposition 2.2. Introducing an auxiliary variable $r=\left(J_{V}^{\varphi^{-1}}\right)^{-1}(\ell-Bu_{n})$, this can be reformulated as a mixed method: find $u_{n}\in U_{n}$ and $r\in V$ such that

[TABLE]

Note that $J_{V}^{\varphi^{-1}}$ is a non-linear map unless $V$ is a Hilbert space. If $V$ is a Hilbert space and we select $\psi(t)=t^{2}/2$, then $J_{V}^{\varphi^{-1}}=R_{V}$ and we recover the framework in [19].

3.2 The Optimal Test Norm

In applications, we are often more interested in minimising the error in the approximation with respect to the norm on $U$ rather than minimising the residual. In other words, we usually want to choose the norm on $V$ in such a way that we ultimately control the error $\|u-u_{n}\|_{U}$. With this in mind, employing the inf-sup condition (3.1a) and continuity of $B$, we deduce that

[TABLE]

where $M$ is the continuity constant of the bilinear form $b$. Hence, in order to control $\|u-u_{n}\|_{U}$, we require $\gamma$ and $M$ as close to one as possible and independent of certain problem specific parameters (the scaling of the diffusion term, for example, in case of the convection-diffusion-reaction equation). The optimal test norm is a concept introduced in the context of DPG methods [37] but unlike the concept of optimal test functions and spaces it can easily be extended to Banach spaces. The optimal test norm is defined as the norm on $V$ such that $M=\gamma=1$ and is given by

[TABLE]

Indeed,

[TABLE]

Conversely, we have the optimal constants $M=\gamma=1$ for any given test norm, if we endow $U$ with the so-called energy norm,

[TABLE]

Ideally, we want to select the norm on $U$ in order to approximate the solution in the desired norm and then work with the optimal test-norm. However, the optimal test-norm is a dual norm which is, in general, not computable. Thus, we have to replace the optimal test-norm with an equivalent norm that is computable. The difficulty here is to obtain equivalence constants, i.e., constants $\gamma$ and $M$, that are independent of problem parameters.

In the context of the convection-diffusion-reaction equation, Broersen and Stevenson work directly with the optimal test norm in [20], whereas the analysis presented in [24] and [22] relies on robust estimates for $\gamma$ and $M$ by looking at the adjoint problem. Extending any of the robust estimates from Hilbert subspaces of $L^{2}\left(\Omega\right)$ to Banach spaces $L^{q}\left(\Omega\right)$, $1<q<\infty$, is highly non-trivial and to the best of our knowledge remains an open problem.

3.3 The Inexact Method

So far, (3.5) does not define a method that can be implemented because the mixed method is still an infinite dimensional problem that relies on the whole test space $V$. Replacing the space $V$ in (3.5) by a finite dimensional $V_{m}\subset V$, however, allows us to approximate the solution to the best-approximation problem by the solution of a finite-dimensional saddle point problem. More precisely, we obtain the following fully discrete mixed problem: find $(r_{m},u_{n})\in V_{m}\times U_{n}$ such that

[TABLE]

where $J_{V}^{\varphi^{-1}}=R_{V}$ if $V$ is a Hilbert space.

Replacing $V$ by a finite-dimensional subspace we obviously lose the best approximation properties. The question of quantifying this ‘loss’ was addressed in [38] for Hilbert spaces by introducing the concept of a Fortin-operator and was later extended to Banach spaces in [27]. For well-posedness of the inexact method we require the existence of a bounded projection operator $\Pi:V\rightarrow V_{m}$ such that

[TABLE]

The norm of $\Pi$ then enters into the a priori estimate and quantifies the ‘loss’ due to the inexactness, i.e., we can obtain an a priori estimate of the form

[TABLE]

Remark 3.1*.*

The constant $C$ can be improved both in the Hilbert setting, where we can obtain $C=\|\Pi\|M/\gamma$, cf., [39], and in the Banach setting, where we would have to introduce geometric constants for the improved estimate, cf., [27].

4 The Convection-Diffusion Equation

The main focus of this article is to apply the abstract approach described in the previous section to the scalar convection-diffusion-reaction equation in a $W^{1,q}\left(\Omega\right)$-$W^{1,q^{\prime}}\left(\Omega\right)$ setting. To this end, consider the following model problem:

[TABLE]

where $\varepsilon$, $\bm{b}:\Omega\rightarrow\mathbb{R}^{d}$ and $c:\Omega\rightarrow\mathbb{R}$ are the (positive) diffusion parameter, convection field and reaction coefficient, respectively, and $f:\Omega\rightarrow\mathbb{R}$ is a given source term. Multiplying by a test function $v\in C_{c}^{\infty}(\Omega)$ and integrating by parts yields the bilinear form

[TABLE]

This allows us to state the following variational problem: find $u\in U:=W^{1,{q}}_{0}\left(\Omega\right)$ such that

[TABLE]

In [28], a proof of an inf-sup condition is given, where the norms on $U$ and $V$ are weighted versions of the standard norms $W^{1,q}\left(\Omega\right)$ and $W^{1,q^{\prime}}\left(\Omega\right)$, respectively. Furthermore, this proof requires certain regularity assumptions on the solution to the Poisson problem and that the convection field $\bm{b}$ and the reaction coefficient $c$ satisfy $c-\nabla\cdot\frac{1}{q}\bm{b}\geq c_{0}>0$. The continuity constant and the inf-sup constant that are established for the bilinear $\mathcal{B}_{\varepsilon}$ in [28] depend on the problem specific parameters. It should be noted that the estimates in [28] can be expected to be sub-optimal since sharper bounds are known for $q=2$.

4.1 Choices for the Test Norm

The choice of the norm for the space $V$ crucially defines the method described in Section 3.3. To this end, we endow $V$ with the following norm

[TABLE]

where $\omega(x)$ is a positive and smooth weighting function. It is well known that for $q=2$ and $\omega(x)\equiv\mathrm{const}\geq 0$ the above choice of norm on $V$ does not yield a robust formulation, i.e., $M/\gamma$ depends on the problem parameters. Here, $M$ denotes the continuity constant of the bilinear form and $\gamma$ its inf-sup constant. The inf-sup constant for $1<q<\infty$ obtained in [28] corresponds to the choice $\omega(x)\equiv 0$ and $\alpha=c_{0}$, where $c_{0}$ is the constant in the Friedrich’s positivity assumption, $c-\frac{1}{q}\nabla\cdot\bm{b}\geq c_{0}$. As mentioned above, the continuity and inf-sup constants are not parameter independent but can be expected to be suboptimal.

One approach to address the robustness issue is to introduce a non-constant weighting function, cf., [24, 19]. For example, we may require $\omega(x)$ to be of magnitude $\mathcal{O}(\varepsilon)$ near the inflow boundary but $\mathcal{O}(1)$ elsewhere. We will show the effect of this by considering a simple one dimensional example on $\Omega=(0,1)$ with the inflow boundary at [math]. For this example $\omega(x)=x+\varepsilon$ has the desired properties; for comparison, we will also consider $\omega(x)\equiv 0$ and $\omega(x)\equiv 1$.

4.2 Weak Boundary Conditions on $r$

Alternatively, the robustness issue can be addressed by changing the boundary conditions on $r$; this has been considered in different ways in [23, 19]. In [22], the boundary conditions on $u$ were modified instead. We will now present the approach in [19] and extend it for $1<q<2$. The idea is to relax the boundary condition on the test space on the inflow part of the boundary. The reasoning behind this is that in the mixed method we essentially approximate the adjoint equation in the test space in order to approximate the residual or the optimal test functions. Under resolved layers at the inflow boundary —which is the outflow boundary for the adjoint equation— then pollute the solution to the primal problem in the inexact method.

Instead of $W^{1,q^{\prime}}_{0}\left(\Omega\right)$, we consider the modified test space $V=W^{1,q^{\prime}}_{0,\Gamma_{+}}(\Omega)$, where $q^{\prime}=q/(q-1)$, and

[TABLE]

Here,

[TABLE]

where $\bm{n}(x)$ denotes the unit outward normal at a point on the boundary $\partial\Omega$. The modified bilinear form (cf., [19]) is given by

[TABLE]

The boundary term is merely the term that is picked up from the integration by parts if $v$ is non-zero on the inflow boundary. Note, however, that the term $\varepsilon\int_{\Gamma_{-}}\frac{\partial u}{\partial\bm{n}}v\mathop{}\!\mathrm{d}s$ is a variational crime if we assume $u\in W^{1,q}_{0}\left(\Omega\right)$. As noted in [19], the correct way of including this term would be to introduce it as an additional unknown on the boundary and the inclusion of $\varepsilon\int_{\Gamma_{-}}\frac{\partial u}{\partial\bm{n}}v\mathop{}\!\mathrm{d}s$ can be viewed as a discrete elimination of this unknown, cf., [19, 20].

4.3 The Inexact Method for the Convection-Diffusion-Reaction Equation

The key step for implementing the inexact method (3.10) for any specific problem is determining the duality mapping $J_{V}^{\varphi^{-1}}$. We choose the weight function $\varphi^{-1}(t)=t^{q^{\prime}-1}$ and compute the duality mapping, similar to Section 2.2; thereby, we get

[TABLE]

For a given right hand side $\ell\in{V}^{\prime}$, we solve the following non-linear system: find $(u_{n},r_{m})\in U_{n}\times V_{m}\subset U\times V$ such that

[TABLE]

One can easily implement (4.8) in, e.g., FEniCS [25, 26] using standard $H^{1}\left(\Omega\right)$-conforming Lagrange finite elements for both $U_{n}$ and $V_{m}$. The spaces $U_{n}$ and $V_{m}$ are chosen over a common mesh. For $U_{n}$ a global polynomial degree $p_{n}$ is chosen and for $V_{m}$ we choose an enriched finite element space with polynomial degree $p_{m}=p_{n}+\varDelta p$, $\varDelta p\geq 1$. To consider the weak boundary conditions on $r$ introduced in the previous section, we simply replace $\mathcal{B}_{\varepsilon}$ with $\tilde{\mathcal{B}}_{\varepsilon}$ and adjust the space $V_{m}$ to only satisfy Dirichlet boundary conditions on the outflow boundary.

4.4 The Limit Case $\varepsilon=0$

Since we are interested in the convection-dominated case, i.e., $\varepsilon\ll\|\bm{b}\|_{L^{\infty}\left(\Omega\right)}$, it makes sense to consider the limit $\varepsilon\rightarrow 0$. Simply setting $\varepsilon=0$ in (4.2) does not yield a well-posed problem unless we only consider boundary conditions on the inflow boundary $\Gamma_{-}$. However, even the ill-posed problem with Dirichlet boundary conditions on the whole boundary $\partial\Omega$ is of interest, since this is essentially the problem that is approximated numerically on coarse meshes when $\varepsilon\ll\|\bm{b}\|_{L^{\infty}\left(\Omega\right)}$, cf., the discussion in [7]. On a discrete level we can think of imposing boundary conditions for the approximation $u_{n}$ on $\Gamma_{+}$ simply as considering an approximation problem in a smaller subspace of $U$.

When $\varepsilon=0$, there are two different weak formulations of the convection-reaction equation depending on whether $\int_{\varOmega}(\bm{b}\cdot\nabla u)v\mathop{}\!\mathrm{d}\bm{x}$ is integrated by parts or not. The two cases differ in the regularity of the trial and test spaces which is reflected in the norms chosen on $U$ and $V$. Thus, in each case the residual is measured in a different norm and (quasi-)best approximation of the analytical solution is achieved in a different space. Both weak formulations can be extended to weak formulations of the convection-diffusion-equation by adding the diffusion term and accounting for different boundary conditions and regularity requirements. If $\varepsilon>0$, we require $u\in W^{1,q}_{0}\left(\Omega\right)$ in both cases and the two weak formulations are formally equivalent. However, the differences observed in the limit case $\varepsilon=0$ can be reflected in the choice of the norm on $W^{1,q}_{0}\left(\Omega\right)$ by choosing the weighting of the terms $\|w\|_{L^{q}\left(\Omega\right)}$, $\|\nabla w\|_{L^{q}\left(\Omega\right)}$ and $\|\nabla\cdot\bm{b}w\|_{L^{q}\left(\Omega\right)}$ accordingly.

We show that for the first choice the exact mixed method is equivalent to the formulation used in [7] and that in the second case a quasi-best approximation in $L^{q}\left(\Omega\right)$ can be computed by determining the optimal test norm. We use this to interpret certain choices of the test norm for the convection-diffusion-reaction equation.

4.4.1 Residual Minimisation in $L^{q}\left(\Omega\right)$

Both weak formulations for the convection-reaction equation are obtained by multiplying the convection-diffusion-reaction equation (4.1) with $\varepsilon=0$ by a smooth test function and integrating over the domain $\Omega$. This immediately yields the first possible choice for the weak formulation: find $u\in W_{0,\Gamma_{-}}^{q}(\bm{b},\varOmega)$ such that

[TABLE]

Here, the bilinear form $\hat{\mathcal{B}}_{0}$ is well defined for $v\in L^{q^{\prime}}\left(\Omega\right)$ and $u$ in the graph space

[TABLE]

endowed with the norm

[TABLE]

The associated residual minimisation problem is

[TABLE]

this was considered in [7]. In general, method (3.10) yields a formulation for solving the residual minimisation problem inexactly. However, since in this case $J^{\varphi}_{V^{\prime}}=J_{q}$, we can avoid the inexactness by directly implementing (3.4); this is exactly the approach considered in [7] where a regularisation is introduced to compute (3.4). The norm on $V=L^{q^{\prime}}\left(\Omega\right)$ corresponds to the choice of the test-norm for the convection-diffusion-reaction equation given in (4.4) with $\omega(x)\equiv 0$ and $\alpha=1$ since this yields the $L^{q^{\prime}}\left(\Omega\right)$-norm if $\varepsilon$ is set to zero. This suggests, that this choice of $\omega(x)$ and $\alpha$ corresponds to the natural extension of this bilinear form to the case $\varepsilon>0$. Moreover, the equivalence of the exact minimum residual method and the approach in [7] for $\varepsilon=0$ together with $\varepsilon$ weighting of the norm of the gradient suggest that for this choice of the test norm, our scheme closely resembles the approach in [7] for $0<\varepsilon\ll 1$.

4.4.2 Best-$L^{q}$ Approximation

Next, we can integrate the convection term by parts in order to obtain the second weak formulation. To this end, we note that

[TABLE]

The boundary term on $\Gamma_{-}$ vanishes since $u=0$ on $\Gamma_{-}$ (for non-zero boundary conditions this term would be absorbed into the right hand side by inserting the boundary condition). For $v$ in the graph space

[TABLE]

the boundary term also vanishes on $\Gamma_{+}$ and we obtain the bilinear form

[TABLE]

which is well-defined for $u\in L^{q}\left(\Omega\right)$. In this case we can compute the optimal test norm

[TABLE]

Choosing the optimal test norm, allows us to obtain the $L^{q}$-best approximation in a given finite dimensional space $U_{n}$. It is easy to see that the optimal test norm is bounded from above up to a constant by the graph norm

[TABLE]

In [40] this formulation of the convection-reaction equation is analysed in detail and an inf-sup condition is established assuming $c-\frac{1}{q}\nabla\cdot\bm{b}\geq c_{0}>0$. This implies that the optimal test norm is up to a constant also bounded from below by the graph norm. In other words, the graph norm is equivalent to the optimal test norm. The graph norm can be obtained from (4.4) by choosing $\omega(x)\equiv 1$, $\alpha=1$ and setting $\varepsilon$ to zero. This suggests that this choice of $\alpha$ and $\omega(x)$ yields the natural extension of this formulation to the convection-diffusion-reaction equation. The dependence of the inf-sup constant in [40] on $c_{0}$ suggests that the equivalence constants can be improved by choosing $\alpha=c_{0}$.

5 Numerical Examples

In this section we consider a range of numerical test cases to illustrate the performance of our proposed numerical scheme. To this end, in Section 5.1 we demonstrate that in the diffusion-dominated regime optimal convergence rates are achieved and moreover that in the convection-dominated regime the convergence rate is as expected. We also study the effect of different choices of the discrete test space $V_{m}$ on the convergence rates. In Section 5.2, we then compare different choices for the test norm on $V$ and the boundary conditions as described in Sections 4.2 and 4.1 for a one-dimensional example. Furthermore, for the same example we study for which choices of $V$, $\|\cdot\|_{V}$ and $V_{m}$ we can observe vanishing over- and undershoots as $q\rightarrow 1$ and whether the method is robust in $\varepsilon$. In Section 5.3 we consider three two-dimensional examples, where for certain selected meshes the overshoots disappear as $q\rightarrow 1$; this behaviour can be predicted by considering $L^{1}$-best approximations of discontinuities.

The four examples we are using in this section are given below. They consist of one simple one-dimensional example and three two-dimensional examples each with solutions containing boundary layers for small $\varepsilon$. The solution to the last example additionally contains an interior layer.

Example 5.1**.**

[TABLE]

The analytical solution is given by

[TABLE]

Example 5.2** (Eriksson-Johnson model problem).**

[TABLE]

The analytical solution is given by

[TABLE]

where

[TABLE]

Example 5.3** (Boundary Layer in the Corner of the Domain).**

[TABLE]

with $\bm{b}=(b_{1},b_{2})^{T}$, $b_{1,2}>0$ and

[TABLE]

The analytical solution is given by

[TABLE]

Since $b_{1,2}>0$, the outflow boundary is defined by the two lines $x=1$ and $y=1$. For small $\varepsilon$ we can observe a boundary layer at the outflow boundary and in particular near the corner $(x,y)=(1,1)$.

Example 5.4** (Interior and Boundary Layer).**

[TABLE]

with $\bm{b}=(2,1)^{T}$.

For this example, a boundary layer develops at $\Gamma_{+}\cap\{y>0.5\}$ and an interior layer along the line $y=0.5x$.

5.1 Convergence Tests

We start with investigating the convergence of the proposed method with weak boundary conditions on the inflow boundary in the space $V_{m}$, i.e., formulation (4.3) with $\tilde{\mathcal{B}}_{\varepsilon}$ from (4.6), and the norm (4.4) with $\omega(x)\equiv 1$ and $\alpha=1$. To this end, we first consider the convergence of the method for Example 5.2 with $\varepsilon=1$. We consider a uniform mesh with element size $h$. The space $U_{n}$ consists of piecewise polynomials of degree $p_{n}$, while space $V_{m}$ consists of piecewise polynomials of degree $p_{m}=p_{n}+\varDelta p$ on the same mesh with $\varDelta p\geq 1$. Figure 1 shows that for $q=1.2$ the error is essentially independent of the choice of $\varDelta p$ and hence there is no benefit to enriching the space of the residual $r_{m}$ beyond $p_{m}=p_{n}+1$. Figure 2 shows that we obtain the expected optimal convergence rates for different polynomial degrees $p_{n}$ under uniform $h$-refinement.

Next, we consider the same example with the same finite element spaces for $\varepsilon=10^{-4}$. Figure 3 (left) shows that again both for $q=1.01$ and $q=1.2$ the error does not improve as $\varDelta p$ is increased beyond $\varDelta p=2$; for $q=1.01$, however, we can observe that the error is larger for $\varDelta p=1$. Note that we did not observe this in the diffusive regime. The plot on the right shows that in the convection-dominated regime we obtain a convergence rate of approximately $\mathcal{O}(h^{\frac{1}{q}})$ as $h$ tends to zero; this is consistent with the approximation error bound of the piecewise linear interpolant of a jump discontinuity.

5.2 Convection-Diffusion in 1D

In this section we consider different choices for $V$, $\|\cdot\|_{V}$ and $V_{m}$ and show how this affects the approximation to the solution $u$ of Example 5.1. To this end, in Section 5.2.1, we investigate how the choice of $\alpha$, $\omega(x)$ and the boundary conditions in $V_{m}$ affect the approximation $u_{n}$. In Section 5.2.2 we consider one specific choice for the test norm and show that the oscillations vanish entirely as $q\rightarrow 1$. Next, in Section 5.2.3, we investigate the choice of the finite element space $V_{m}$ more closely. We have seen that increasing $\varDelta p$ has no significant effect on the convergence rates; however, we will see that this does indeed affect certain qualitative properties of the solution. Finally, we demonstrate that the method is robust in $\varepsilon$ in Section 5.2.4

5.2.1 Comparing the Choices for the Test Norm and the Boundary Conditions

We start by comparing different versions of the non-linear Petrov-Galerkin method. In order to ensure that the space $V_{m}$ is sufficiently large we use piecewise polynomials of degree $p_{m}=10$ as the basis, while $U_{n}$ consists of piecewise linear polynomials. We split the interval $(0,1)$ uniformly into $8$ elements. We consider three different choices for the weighting function $\omega(x)$ as mentioned in Section 4.1, namely $\omega(x)=x+\varepsilon$, $\omega(x)\equiv 1$ and $\omega(x)\equiv 0$. We combine this with two different choices of boundary conditions on $r$, i.e., zero Dirichlet boundary conditions on the whole boundary and weak boundary conditions on the inflow boundary as described in Section 4.2. This creates six test cases; we first consider these six test cases for $\varepsilon=10^{-3}$ with $\alpha=1$ and $q=2,1.01$. Secondly, we consider all six test cases with $\varepsilon=10^{-6}$, $\alpha=0,1$ and $q=2,1.01$ The solution $u_{n}$ for each of the cases is shown in Figure 4.

We can see in Figures 4(a), 4(c) and 4(e) that for $q=2$ and Dirichlet boundary conditions on the whole boundary for $r$, the approximation $u_{n}$ only resembles the analytical solution if $\omega(x)=x+\varepsilon$. Introducing weak boundary conditions on $r$ on the inflow boundary resolves this issue for $\omega(x)\equiv 1$ but the approximation with $\omega(x)\equiv 0$ only shows improvement if $\alpha=0$. These observations are the same for $\varepsilon=10^{-3}$ and $\varepsilon=10^{-6}$. This indicates that we can circumvent constructing a non-constant function $\omega(x)$ – which can be challenging for complicated geometries – by introducing weak boundary conditions on $r$ on the inflow boundary. These observations are consistent with the results for a very similar problem studied in [19]. Furthermore, note that $\omega(x)=0$ and $\alpha=1$ most closely resembles the method introduced in [7]. In [7] it is demonstrated for a different (two-dimensional) example that the approximation for $q=2$ can be very inaccurate.

If we now consider $q=1.01$, cf., Figures 4(b), 4(d) and 4(f), we observe a different behaviour of the approximation. If $\varepsilon=10^{-3}$, the approximations are all very similar and much closer to the analytical solution with close to no undershoots; we will investigate the phenomenon of vanishing undershoots more closely both in one dimension and two dimensions later on. If $\varepsilon=10^{-6}$, the approximation with Dirichlet boundary conditions on the whole boundary and $\omega(x)\equiv 0$ or $\omega(x)\equiv 1$ again leads to a very poor approximation to the analytical solution. This again points to robustness issues that for $q=1.01$ are only visible for much smaller $\varepsilon$ than for $q=2$. An interesting observation is that for $\alpha=1$, the combination of weak inflow boundary conditions on $r$ and $\omega(x)\equiv 0$ is not distinguishable from other choices of $\omega(x)$ for $q=1.01$ whereas, for $q=2$ it is significantly different. The improvement of the approximation as $q\rightarrow 1$ for $\omega(x)\equiv 0$ and $\alpha=1$ resembles the results presented in [7].

In conclusion, choosing the boundary conditions as described in Section 4.2 allows us to avoid constructing a non-constant weighting function $\omega(x)$. Close to $q=1$, the results are very similar for $\omega(x)\equiv 1$ and $\omega(x)\equiv 0$, whereas close to $q=2$, $\omega(x)\equiv 1$ clearly yields a significantly better approximation of the solution.

The choice of $\alpha$ does not seem to have a big impact on the approximation except in the case $\omega(x)\equiv 0$ and $q=2$, where choosing $\alpha=0$ improves the approximation compared to $\alpha=1$. If $\omega(x)\equiv 1$, then both for $\alpha=1$ and $\alpha=0$, the dominant term in the norm on $V$ is essentially $\|\nabla u\|_{L^{q^{\prime}}\left(\Omega\right)}$ since $\|\bm{b}\cdot\nabla v\|_{L^{q^{\prime}}\left(\Omega\right)}$ and $\|\nabla v\|_{L^{q^{\prime}}\left(\Omega\right)}$ are the same up to a constant in one dimension. If $\omega(x)\equiv 0$, then $\alpha=0$ implies that $V$ is endowed with the $W^{1,q^{\prime}}_{0}\left(\Omega\right)$-norm , whereas for $\alpha=1$ the $L^{q^{\prime}}\left(\Omega\right)$-term is the dominant term in the norm on $V$ due to the $\varepsilon$-weighting of the gradient. Since $\alpha=0$ yields a better approximation in this case, this suggests that for $q=2$ it is favourable to choose a stronger norm on $V$ that is either dominated by the norm of the gradient or the streamline term $\|\bm{b}\cdot\nabla v\|_{L^{q^{\prime}}\left(\Omega\right)}$.

5.2.2 Vanishing oscillations as $q\rightarrow 1$

The examples in Figure 4 already illustrate that the undershoot in the approximate solution nearly vanishes for $q=1.01$. We will now investigate in more detail how the undershoot depends on $q$ if we choose $\varepsilon=10^{-5}$, $\omega(x)\equiv 1$, $\alpha=1$ and impose weak inflow boundary conditions on $r_{m}$. To this end, we choose a large space $V_{m}$ with polynomial degree $p_{m}=10$ and piecewise linear polynomials for $U_{n}$. We split the interval into $8$ elements and compute the approximate solutions $(r_{m},u_{n})$ for several choices of $q$. Figure 5 shows $u_{n}$ on the left and $r_{m}$ on the right. A close-up of the undershoot and a plot $\min(u_{n})$ vs. $q$ is shown in the centre of the figure. We can see that the undershoot decreases monotonically as $q$ approaches $1$. In contrast to this, we can see that the oscillations in the residual $r_{m}$ increase as $q^{\prime}\rightarrow\infty$ which suggests that a large space $V_{m}$ may be necessary. This is the next aspect of the method we will investigate.

5.2.3 The Choice of the Space $V_{m}$

In [17], it has been shown that the inexact method is well-posed if a Fortin projector, cf., (3.11) exists. It is easy to see that $\dim(V_{m})\geq\dim(U_{n})$ is a necessary condition. Finding a sufficient condition or in other words finding a compatible pair $(U_{n},V_{m})$ is highly non-trivial. In [40], certain special cases for the convection-reaction equation are considered. The observations in Section 5.1 suggest that for $\varDelta p\geq 2$, the Fortin condition is typically satisfied. We now investigate how the choice of $V_{m}$ affects the undershoot for $q$ close to $1$. To this end, we consider two different strategies of enlarging $V_{m}$: global $p$-enrichment and uniform $h$-refinement. We again consider $\varepsilon=10^{-5}$, $\omega(x)\equiv 1$, $\alpha=1$ and impose weak inflow boundary conditions on $r_{m}$.

p-enrichment: In terms of implementation, the simplest choice for the space $V_{m}$ such that $\mathrm{dim}(V_{m})>\mathrm{dim}(U_{n})$ is a finite element space over the same grid with the polynomial degree increased globally by some integer $\varDelta p$. Figure 6 shows how the overshoot reduces as we increase $\varDelta p$. The difference in $r_{m}$ for varying $p_{m}$ is barely visible even though we have used a finer mesh for plotting in order to capture the behaviour of the higher modes and yet it has a significant effect on the approximation $u_{n}$.

h-refinement: This time both $U_{n}$ and $V_{m}$ consist of piecewise linear polynomials. To ensure that $\mathrm{dim}(V_{m})>\mathrm{dim}(U_{n})$, we choose a refinement of the underlying mesh of the space $U_{n}$ to construct the space $V_{m}$. This is the only numerical experiment in this article that cannot be implemented in FEniCS; instead, e.g., the C++ library Hermes2D [41] can be used. Figure 7 shows how the overshoot reduces as we refine the mesh for $V_{m}$;

the difference in $r_{m}$ is again barely visible despite the significant effect on the approximation $u_{n}$.

5.2.4 Robustness in $\varepsilon$

As a final experiment in one dimension, we study how $\varepsilon$ affects the undershoot in the approximation. We continue with the same setting as in the previous section, but keep $V_{m}$ fixed as the space of piecewise polynomials of degree $p_{m}=10$ on the same mesh as used for $U_{n}$ and vary $\varepsilon$ instead. Figure 8 shows that, although the undershoot is not the same for all $\varepsilon$ leading to an under resolved layer, the undershoot seems to be approximately the same for all $\varepsilon\leq\varepsilon_{0}$ for some $\varepsilon_{0}$. This can be traced back to the inflow boundary conditions on $r_{m}$ — which are scaled with $\varepsilon$ — if we look at $r_{m}$ on the right in Figure 8. The method therefore seems to be robust in $\varepsilon$.

5.3 Vanishing Oscillations in Two Dimensions

In this section we explore whether the under- and overshoots in the approximation still vanish as $q\rightarrow 1$ if we apply the method to two-dimensional problems. To this end, we consider the Examples 5.2, 5.3 and 5.4 on different meshes for $\varepsilon=10^{-6}$. We consider $\omega(x)\equiv 1$, $\alpha=1$ and impose weak boundary conditions on $r_{m}$ as described in Section 4.2. In one case we will also consider $\alpha=0$ to compare the approximations for both choices of $\alpha$. We will see that the overshoots disappear on certain meshes, but remain present on others. We will furthermore demonstrate that the approximation qualitatively behaves like the $L^{q}(\Omega)$-best approximation of the analytical solution and thus we can apply the observations in [5] to (a) predict on which meshes the overshoot will disappear as $q\rightarrow 1$ and (b) design meshes that have this property in specific situations.

5.3.1 The Eriksson-Johnson Model Problem

We consider the Eriksson-Johnson Model Problem (Example 5.2) on the four different meshes depicted in Figure 9. Note that these are the same meshes that are investigated in [5].

First, we compare the approximations on Mesh 2 and Mesh 3 for $q=2$ and $q=1.01$ with the interpolant of the analytical solution in $U_{n}$. Figure 10 shows the approximations on Mesh 2; for $q=2$ we clearly observe overshoots near the boundary layer (center) compared to the interpolant of the analytical solution (left). On the other hand, for $q=1.01$ we can see that the approximation approaches the interpolant of the analytical solution. Figure 11 shows the approximations on Mesh 3; in this case we can see that the overshoot is nearly the same for $q=2$ and $q=1.01$. From [5] we note that this is the qualitative behaviour that the $L^{q}$-best approximation exhibits.

With this in mind, we have numerically computed the $L^{q}(\Omega)$-best approximation of the analytical solution for several choices of $q$ and have also computed the finite element approximations with $\omega(x)=1$ and both, $\alpha=1$ and $\alpha=0$, with weak boundary conditions imposed on $r_{m}$ on the inflow boundary for the same values of $q$. Figure 12 shows the maximal error between the interpolant of the analytical solution and the approximation for the $L^{q}(\Omega)$-best approximation and the error between the interpolant and the numerical approximations for $\alpha=1$ and $\alpha=0$. We can clearly see that the overshoot observed for our proposed scheme is very similar to the overshoot in the $L^{q}(\Omega)$-best approximation. Therefore, understanding for which meshes discontinuities/layers are captured as sharply as the grid allows by the $L^{q}(\Omega)$-best approximation can serve as an indicator regarding whether the over- and undershoots dissappear as $q\rightarrow 1$. This is the subject of [5], where the $L^{q}$-best approximation of discontinuous functions is analysed in detail in certain cases.

5.3.2 Boundary Layer in a Corner of the Domain

In this section we consider Example 5.3 with $\bm{b}=(2,1)^{T}$. We will see that on Mesh 2 in this case the overshoot does not disappear as $q\rightarrow 1$ for our method. This is due to the layer not only appearing along an edge of the unit square but also near the corner $(x,y)=(1,1)$. With this in mind, we suggest an alternative mesh that is only a slight modification of Mesh 2 for which the overshoot disappears. Figure 13 shows the two versions of Mesh 2 that we have used and the corresponding approximations for $q=2$ and $q=1.01$ along the line through the interior nodes closest to the boundary $y=1$. We can see that on Mesh 2, the overshoot near the corner does not disappear as $q\rightarrow 1$, whereas it does reduce significantly away from the corner. For the modified version of Mesh 2, the overshoot disappears everywhere as $q\rightarrow 1$.

To understand the improvement when using this modified mesh, we focus on the interior node closest to the corner $(x,y)=(1,1)$. Figure 14 shows all elements connected to this node. The elements with either a node or an edge on the boundary are marked in green; the elements separated from the boundary are blue. From [5] we can infer that the $L^{1}$-best approximation does not exhibit overshoots if the volume of the green area is smaller or equal to the volume of the blue area. This is obviously violated on Mesh 2 with the green area being three times as large as the blue area. The modified version of the mesh is designed to satisfy this condition by changing the mesh as indicated in Figure 14.

5.3.3 Example with an Interior Layer

In this section we consider Example 5.4 and demonstrate that a mesh can be constructed such that the over- and undershoots both near the boundary layer and the interior layer can be eliminated as $q\rightarrow 1$. To this end, Figure 15 shows the approximation to Example 5.4 for $q=2$ and $q=1.1$ on two different meshes. The first mesh is designed such that the overshoot at the boundary layer disappears as $q\rightarrow 1$; this is already clearly visible for $q=1.1$. Along the interior layer, the overshoot also reduces significantly, but some over- and undershoots are still present for $q=1.1$. The second mesh is additionally designed to align with the interior layer in such a way that we can expect the over- and undershoots to disappear as $q\rightarrow 1$. We can see, that this mesh even nearly eliminates the overshoot along the interior layer for $q=2$ and for $q=1.1$ all over- and undershoots along both layers have essentially vanished.

Figure 16 shows the approximations for $q=2$ and $q=1.1$ if $\bm{b}$ is not perfectly aligned with the mesh, i.e., $\bm{b}=(2,1.2)^{T}$ and $\bm{b}=(2,1.06)^{T}$. The difference between these two choices for $\bm{b}$ is that in the latter case the layer is still contained in between the two lines parallel to the line $(x,y)=t(2,1)$, whereas in the former case it is not. In both cases, we clearly observe overshoots along the interior layer for $q=2$ which are reduced for $q=1.1$. Both, for $q=2$ and $q=1.1$ the overshoot reduces slowly the closer $\bm{b}$ is to $(2,1)^{T}$. For $\bm{b}=(2,1.06)^{T}$ and $q=1.1$ the overshoot is barely visible in Figure 16, but comparing the maximum and minimum values of the approximations reveals that the over- and undershoots are similar to the case $\bm{b}=(2,1.2)$ in magnitude.

6 Conclusions and Future Directions

In this article, we generalized the framework employed for DPG methods to more general Banach spaces. This framework in principle allows us to choose the norm in which the solution is approximated, provided a suitable well-posed variational formulation with robust constants and a test-norm that is computable is given. Another challenge in designing such a method is that a Fortin-condition must be satisfied for the discrete spaces $U_{n}$ and $V_{m}$. The recent work [42] describes a double adaptivity approach for the DPG method to circumvent this challenge in the context of DPG methods. Due to the close relationship between DPG methods and our approach this has potential to be generalized to the Banach space framework.

Furthermore, we have demonstrated how to design a finite element method for the convection-diffusion-reaction problem that in the convection-dominated case yields solutions that qualitatively behave like the $L^{q}(\Omega)$-best approximation of the analytical solution. This means that on meshes where the $L^{1}(\Omega)$-best approximation does not contain any over- and undershoots, these can also be avoided by taking $q\rightarrow 1$ in our proposed method. Indeed, the final two examples presented in this article demonstrate that it is possible to adjust meshes to specific situations in order to eliminate over- and undershoots.

An important implication of the connection to the $L^{q}\left(\Omega\right)$-best approximation is that we can use the insights from [5] to design suitable meshes. Note that even in one dimension the over- and undershoots will not disappear for certain non-uniform meshes. Roughly speaking, meshes where the elements near the discontinuity or layer present in the analytical solution are smaller than at a distance from those features are more favourable. Additionally, it is desirable that all interior nodes that are closest to such a feature are aligned roughly in parallel with the feature; for further details, we refer to [5]. This is similar to the observations used to design the modified version of the mesh in the final two examples. This suggests that if the mesh is refined in a suitable way near the layer, we can eliminate over- and undershoots. This observation could be used to design a mesh refinement strategy that modifies the mesh in such a way that the $L^{1}(\Omega)$-best approximation of the layer or discontinuity does not contain any over- or undershoots. Note, however, that even though we modified the mesh in the final example in a way that resembles refinement towards the boundary layer, the boundary layer clearly remains under resolved even on the modified mesh. It is therefore not necessary to fully resolve the layer to eliminate the oscillations in this case.

Another challenge of our proposed method is computational feasibility. Possible strategies to addressing this include reducing the degrees of freedom by choosing the test space $V_{m}$ adaptively, introducing broken test spaces $V_{m}$, similar to DPG methods, to allow for more effective parallelization, and developing a more efficient non-linear solver. Despite the challenges outlined above, the numerical results for our proposed method show great potential for the application to non-linear problems whose analytical solutions contain sharp layers or discontinuities. It is well-known that higher-order monotonicity preserving methods can be expected to be non-linear; thus, the non-linear nature of our approach by no means diminishes its potential.

Acknowledgements

The authors would like to thank Gabriel Barrenechea, Jesse Chan, Leszek Demkowicz, Jay Gopalakrishnan, Ignacio Muga and Martin Stynes for fruitful discussions on the topic of this article.