The Discrete-Dual Minimal-Residual Method (DDMRes) for Weak Advection-Reaction Problems in Banach spaces

We propose and analyse a minimal-residual method in discrete dual norms for approximating the solution of the advection-reaction equation in a weak Banach-space setting. The weak formulation allows for the direct approximation of solutions in the Lebesgue $L^p$-space, $1<p<\infty$. The greater generality of this weak setting is natural when dealing with rough data and highly irregular solutions, and when enhanced qualitative features of the approximations are needed. We first present a rigorous analysis of the well-posedness of the underlying continuous weak formulation, under natural assumptions on the advection-reaction coefficients. The main contribution is the study of several discrete subspace pairs guaranteeing the discrete stability of the method and quasi-optimality in $L^p$, and providing numerical illustrations of these findings, including the elimination of Gibbs phenomena, computation of optimal test spaces, and application to 2-D advection.


Introduction
Residual minimization encapsulates the idea that an approximation to the solution u ∈ U of an (infinite-dimensional) operator equation Bu = f , can be found by minimizing the norm of the residual f − Bw n amongst all w n in some finitedimensional subspace U n ⊂ U. This powerful idea provides a stable and convergent discretization method under quite general assumptions, i.e., when B : U → V * is any linear continuous bijection from Banach space U onto the dual V * of a Banach space V, f ∈ V * , and dist(u, U n ) → 0 as n → ∞; see, e.g., Guermond [29,Section 2] for details. Note that this applies to well-posed weak formulations of linear partial differential equations (PDEs), in which case B is induced by the underlying bilinear form (i.e., Bw, v = b(w, v) , ∀w ∈ U, ∀v ∈ V). As such, residual minimization is essentially an ideal methodology for non-coercive and/or nonsymmetric problems.
However, for many weak formulations of PDEs, V * is a negative space (such as H −m (Ω), or more generally W −m,p (Ω), which is the dual of the Sobolev space W m,q 0 (Ω), where 1 < p < ∞, p −1 + q −1 = 1, m = 1, 2, . . . , or the dual of a graph space). In that case, this requires the minimization of the residual in the non-computable dual norm · V * . To make this tractable, one can instead minimize in a discrete dual norm. In other words, one aims to find an approximation u n ∈ U n such that f − Bu n (Vm) * is minimal, where V m is some finite-dimensional subspace of V. We refer to this discretization method as residual minimization in discrete dual norms, or simply as the DDMRes method (Discrete-Dual Minimal-Residual method).
In this paper, we consider the DDMRes method when applied to a canonical linear first-order PDE in weak Banach-space settings. In particular, we consider the advection-reaction operator u → β ·∇u+µ u, with β : Ω → R d and µ : Ω → R given advection-reaction coefficients, in a functional setting for which the solution space U is L p (Ω), 1 < p < ∞, and V is a suitable Banach graph space (see Section 3 for details). This weak setting allows for the direct approximation of irregular solutions, while the greater generality of Banach spaces (over more common Hilbert spaces) is useful for example in the extension to nonlinear hyperbolic PDEs [30], 1 as well as in approximating solutions with discontinuities (allowing the elimination of Gibbs phenomena; see further details below).
It has recently become clear that many methods are equivalent to DDMRes, the most well-known being the discontinuous Petrov-Galerkin (DPG) method (for which B corresponds to a hybrid formulation of the underlying problem, so that V is a broken Sobolev-type space), see Demkowicz and Gopalakrishnan [20], and the Petrov-Galerkin method with projected optimal test spaces, see Dahmen et al. [15]. While these methods require U and V to be Hilbert spaces, in more general Banach spaces the DDMRes method is equivalent to certain (inexact) nonlinear Petrov-Galerkin methods, or equivalently, mixed methods with monotone nonlinearity, where the nonlinearity originates from the nonlinear duality map J V : V → V * ; see Muga & Van der Zee [32] for details, including a schematic overview of connections to other methods.
The numerical analysis of the DDMRes method has been carried out abstractly by Gopalakrishnan & Qiu [27] in Hilbert spaces (see also [15,Section 3]), and by Muga & Van der Zee [32] in smooth Banach spaces. A key requirement in these analyses is the Fortin compatibility condition on the family of discrete subspace pairs (U n , V m ) under consideration, which, once established, implies stability and quasi-optimal convergence of the method. In some sense, the Fortin condition is rather mild, since for a given U n , there is the expectation that it will be satisfied for a sufficiently large V m (thereby making the discrete dual norm · (Vm) * sufficiently close to · V * ). Of course, whether this can be established depends crucially on the operator B, therefore also on the particular weak formulation of the PDE that is being studied.
The main contribution of this paper consists in the study of several elementary discrete subspace pairs (U n , V m ) for the DDMRes method for weak advectionreaction, including proofs of Fortin compatibility in the above-mentioned Banachspace setting. It thereby provides the first application and corresponding analysis of DDMRes in genuine (non-Hilbert) Banach spaces. In particular, for the given compatible pairs, DDMRes is thus a quasi-optimal method providing a near-best approximation in L p (Ω). Note that our results do not cover DPG-type hybrid weak formulations (with a broken graph space V), so that our discrete spaces V m are globally conforming. Broken Banach-space settings will be treated in forthcoming work.
We now briefly discuss some details of our results. To be able to carry out the analysis, our results focus on discrete subspace pairs (U n , V m ), where U n is a lowest-order finite element space on mesh T n in certain specialized settings. We first consider continuous linear finite elements in combination with continuous finite elements of degree k, i.e., U n = P 1 cont (T n ) and V m = P k cont (T n ). The Fortin condition holds when k ≥ 2, assuming, e.g., incompressible pure advection (div β = µ = 0) in a one-dimensional setting. Interestingly, we demonstrate that the notorious Gibbs phenomenon of spurious numerical over-and undershoots, commonly encountered while approximating discontinuous solutions with continuous approximations, can be eliminated with the DDMRes method upon p → 1 + (see Section 5.1), which is in agreement with previous findings on L 1 -methods [29,31].
We then consider U n = P 0 (T n ), that is, discontinuous piecewise-constant approximations on arbitrary partitionings of the domain Ω in R d , d ≥ 1. It turns out that it is possible to define an optimal test space S n := B − * U n and subsequently prove Fortin's condition for any V m ⊇ S n . This result essentially hinges on the fact that U n is invariant under the L p duality map (see proof of Proposition 5.3). Since the optimal test space is however not explicit, it requires in general the computation of an explicit basis (see Section 5.2). Such computations may not be feasible in practise, and in those cases, as an alternative, one could resort to sufficiently-rich V m , e.g., continuous linear finite elements on a sufficiently-refined submesh of the original mesh (cf. [7]).
Interestingly however, under certain special, yet nontrivial situations, the optimal test space S n happens to coincide with a convenient finite element space. For example, in 2-D in the incompressible pure advection case with β piecewise constant on some partition, if T n is a triangular mesh of Ω (compatible with the partition) and all triangles are flow-aligned, then we prove that S n = P 1 conf (T n ), where P 1 conf (T n ) refers to the space of piecewise-linear functions that are conforming with respect to the graph space V. Numerical experiments in 2-D indeed confirm in this case the quasi-optimality of DDMRes (see Section 5.3).
In recent years, several similar methods for weak advection-reaction have appeared, all of which were in Hilbert-space settings (i.e., the solution space is L 2 (Ω)) and use a broken weak formulation. Indeed, these include some of the initial DPG methods [18,19,8], which were proposed before the importance of Fortin's condition was clarified. Recently however, Broersen, Dahmen & Stevenson [7] studied a higher-order pair using standard finite-element spaces for the DPG method of weak advection-reaction. Under mild conditions on β, they proved Fortin's condition when U n consists of piecewise polynomials of degree k, and V m consists of piecewise polynomials of higher degree over a sufficiently-deep refinement of the trial mesh. The extension of their proof, based on approximating the optimal test space, to any Banach-space setting seems nontrivial, since currently the concept of an optimal test space is in general absent in DDMRes in Banach spaces (cf. [32]), exceptions notwithstanding (such as the lowest-order piecewise-constant case discussed above).
Let us finally point out that methods for weak advection-reaction are quite distinct from methods for strong advection-reaction (which has its residual in L p (Ω) and a priori demands more regularity on its solution). Indeed, there is a plethora of methods in the strong case; see, e.g., Ern & Guermond [22,Chapter 5] and Guermond [29], all of which typically exhibit suboptimal convergence behaviour when measured in L p (Ω). In the context of strong advection-reaction the results by Guermond [28] are noteworthy, who proved the Fortin condition for several pairs, consisting of a low-order finite element space and its enrichment with bubbles. These results however do not apply to weak advection-reaction. Similarly for the stability result by Chan, Evans & Qiu [13].
The remainder of this paper is arranged as follows. In Section 2, we first present preliminaries for the advection-reaction equation, allowing us to recall in Section 3 the specifics of the well-posed Banach-space setting (cf. Cantin [9]). In particular, we provide a self-contained proof of the continuous inf-sup conditions using various properties of the L p duality map. Then, in Section 4, we consider the discrete problem corresponding to the DDMRes method in the equivalent form of the monotone mixed method, and establish stability and quasi-optimality of the method, provided the Fortin condition holds. In Section 5, we consider particular discrete subspace pairs (U n , V m ). This section contains several proofs of Fortin conditions, as well as some illustrative numerical examples pertaining to the Gibbs phenomena (Section 5.1), optimal test space basis (Section 5.2), and quasi-optimal convergence for 2-D advection (Section 5.3).

Advection-reaction preliminaries
For any dimension d ≥ 1, let Ω ⊂ R d be an open bounded domain, with Lipschitz boundary ∂Ω oriented by a unit outward normal vector n. Let β ∈ L ∞ (Ω) be an advection-field such that div β ∈ L ∞ (Ω) and let µ ∈ L ∞ (Ω) be a (spacedependent) reaction coefficient. The advection-field splits the boundary ∂Ω into an inflow, outflow and characteristic part, which for continuous β corresponds to respectively; see [7,Section 2] for the definition of the parts in the more general case β, div β ∈ L ∞ (Ω) (which is based on the integration-by-parts formula (2.7)).
Given a possibly rough source f • and inflow data g, the advection-reaction model is: Before we give a weak formulation for this model and discuss it's well-posedness, we first introduce relevant assumptions and function spaces. We have in mind a weak setting where u ∈ L p (Ω), for any p in (1, ∞). Therefore, throughout this section, let 1 < p < ∞ and let q ∈ (1, ∞) denote the conjugate exponent of p, satisfying the relation p −1 + q −1 = 1.
The following assumptions are natural extensions of the classical ones in the Hilbert case.
Assumption 2.1 (Friedrich's positivity assumption) There exists a constant µ 0 > 0 for which: Assumption 2.2 (Well-separated in-and outflow) The in-and outflow boundaries are well-separated, i.e., ∂Ω − ∩ ∂Ω + = ∅ and, by partition of unity, there exists a function For brevity we use the following notation for norms and duality pairings: endowed with the norm The "adjoint" norm is defined by These norms are equivalent, which can be shown by means of the identity div(βw) = div(β)w + β · ∇w .
(2.6) Remark 2.4 (Graph-spaces and traces) As a consequence of Assumption 2.2, traces on the space W ρ (β; Ω) are well-defined as functions in the space and, moreover, for all w ∈ W ρ (β; Ω) and all v ∈ W σ (β; Ω), the following integration-by-parts formula holds: The proof of these results is a straightforward extension of the Hilbert-space case given in, e.g., Di Pietro & Ern [ . We can thus define the following two closed subspaces, which are relevant for prescribing boundary conditions at ∂Ω + or ∂Ω − :

Remark 2.5 (Non-separated in-and outflow)
The requirement of separated in-and outflow can be removed, but different trace operators have to be introduced [26].
Remark 2.9 (Weaker statements) Under a weaker condition than Assumption 2.1, Lemma 2.8 can be generalized to the following situations: Indeed, it is enough to verify the existence of a constant µ * 0 > 0 and a Lipschitz continuous function ζ(x) such that: These statements can be inferred from the recent work of Cantin [9]. Notice that if Assumption 2.1 is satisfied, then (2.15) holds with ζ(x) ≡ 0 and µ * 0 = µ 0 .

A weak setting for advection-reaction
The weak setting for the advection-reaction problem (2.1) considers a trial space , and a test space V := W q 0,+ (β; Ω) endowed with the norm ||| · ||| q,β (see (2.5)). The weak-formulation reads as follows: and the right-hand side f is related to the original PDE data (f • , g) via: where f • is given in (for example) L p (Ω) and g is given in L p (|β · n|; ∂Ω). More rough f • is allowed as long as f ∈ [W q 0,+ (β; Ω)] * .
The following result states the well-posedness of problem (3.1). Although this result can be inferred from the recent result by Cantin [9], we provide a slightly alternative proof based on establishing the adjoint inf-sup conditions using properties of the L p duality map. For a classical proof of well-posedness in a similar Banach-space setting, we refer to Beirão Da Veiga [4,5] (cf. [3] and [17, Chapter XXI]). (ii) In the case that Assumption 2.1 holds true, we have the following a priori bound: and µ 0 > 0 being the constant in Assumption 2.1.
(ii ) On the other hand, in the case where Assumption 2.6 holds true, we also have the a priori bound Proof See Section A.1.

The general discrete problem
We now consider the approximate solution of (3.1) given by the DDMRes method, i.e., given finite-dimensional subspaces U n ⊂ U = L p (Ω) and V m ⊂ V = W q 0,+ (β; Ω), we aim to find u n ∈ U n such that: where the discrete dual norm is given by As proven in [32,Theorem 4.A], the minimization problem (4.1) is equivalent to following monotone mixed method: In (4.2a), J V : V → V * denotes the (monotone and nonlinear) duality map of V = W q 0,+ (β; Ω) defined by the action: is the duality map of L q (Ω). The solution u n ∈ U n of (4.2) is exactly the residual minimizer of (4.1), while r m ∈ V m is a representative of the discrete residual, i.e., The well-posedness of the discrete method (4.2) relies on the well-posedness of the continuous problem (3.1) (see Theorem 3.A), together with the following Fortin assumption.
where I : V → V is the identity map in V. For simplicity, we write Π instead of Π n,m . (i) There exists a unique solution (r m , u n ) to (4.2), which satisfies the a priori bounds: (ii) Moreover, we have the a priori error estimate:  A implies optimal convergence rates in L p (Ω) for finite element subspaces U n ≡ U h , provided C Π is uniformly bounded. For example, on a sequence of approximation subspaces {P k (T h )} h>0 of piecewise polynomials of fixed degree k on quasi-uniform shape-regular meshes T h with mesh-size parameter h, well-known best-approximation estimates (see, e.g., [6,Section 4.4], [22, Section 1.5] and [24]) imply where | · | W s,p (Ω) denotes a standard semi-norm of W s,p (Ω) (e.g., of Sobolev-Slobodeckij type). For a relevant regularity result in W 1,ρ (Ω), with ρ ≥ 2, see Girault & Tartar [25] (see also [33]).

Applications
In this section we apply the general discrete method (4.2) to particular choices of discrete subspace pairs (U n , V m ) involving low-order finite-element spaces. For simplicity, throughout this section Ω ⊂ R d will be a polyhedral domain and T n will denote a finite partition of Ω, i.e., T n = {T } consists of a finite number of non-overlapping elements T for which Ω = T ∈Tn T . be 2δ 0 , where δ 0 is the Dirac delta at x = 0, i.e., Notice that the exact solution of (5.1) corresponds to the sign of x: We endow V = W q 0,+ (β; Ω) = W 1,q 0,{1} (Ω) with the norm · V = (·) q , which simplifies the duality map (4.3) to a normalized q-Laplace operator : Moreover, in this setting, it is not difficult to show that residual minimization in [W 1,q 0,{1} (Ω)] * now coincides with finding the best L p -approximation to sign(x). The Gibbs phenomena for best L p -approximation was studied analytically by Saff & Tashev [34], when using continuous piecewise-linear approximations on n equalsized elements. They clarified that the overshoot next to a discontinuity remains as n → ∞ whenever p > 1, however remarkably, the overshoot tends to zero as p → 1 + .
To illustrate these findings, we plot in Figure 1 the best L 2 -approximation using continuous piecewise-linears for various mesh-size parameter h. Clearly, the overshoots remain present, signifying the Gibbs phenomenon. Next, in Figure 2 we plot the solution to ideal residual minimization (i.e., in the so-called ideal case where V m = V) on a fixed mesh consisting of nine elements, for different values of p > 1. 2 In Figure 2, we also plot the corresponding ideal residual r (x) as defined by the mixed formulation (4.2) in the case where V m = V. It can be shown that in this ideal 1-D situation r = u n − u 2−p p |u n − u| p−1 sign(u n − u). The plots in Figure 2 clearly illustrate the elimination of the Gibbs phenomenon as p → 1 + . We note that the elimination of the Gibbs phenomena was also observed for residual minimization of the strong form of the advection-reaction residual in L 1 (Ω) [29], the explanation of which remains somewhat elusive.
Next, consider the DDMRes method (4.2) with the discrete space pair (U n , V m ) defined by: where k is the polynomial degree of the test space, and T n is any partition of the interval (−1, 1).
Proof See Section A.2.
Remark 5.2 (Solution of nonlinear system) The DDMRes method leads to a discrete (nonlinear) q-Laplace mixed system, which, for p moderately close to 2, can be solved with, e.g., Newton's or Picard's method. For p close to 1 or much larger than 2 the nonlinear problem becomes more tedious to solve, and we have resorted to continuation techniques (with respect to p) or a descent method for the equivalent constrained-minimization formulation. Figure 3 diplays numerical results obtained using the DDMRes method (4.2) with the above discrete spaces and for p = 1.01 (hence q = 101). We plot u n and r m for various test-space degree k ≥ 2. While the method is stable for any k ≥ 2 (owing to Proposition 5.1), there is no reason for the DDMRes method to directly inherit any qualitative feature of the exact residual minimization (i.e., V m = V). Indeed, overshoot is present for small k. However, results are qualitatively converging once the test space V m starts resolving the ideal r. This is expected by [32,Proposition 4.2], which states that the ideal u n is obtained if the ideal r happens to be in V m . The lines indicated by "ideal" in Figure 3 correspond to the case that r is fully resolved.
Interestingly, the results seem to indicate that different values of k in each element would be needed for efficiently addressing Gibbs phenomena. This is reminiscent of the idea of adaptive stabilization [14], in which for a given U n a sufficiently large test-space V m is found in an adaptive manner so as to achieve stability.

5.2
The pair P 0 (T n ) -B − * (P 0 (T n )): An optimal compatible pair In the remainder of Section 5, we consider for U n piecewise-constant functions on mesh-partitions T n of Ω ⊂ R d , i.e., U n ⊆ P 0 (T n ) := w n ∈ L ∞ (Ω) : w n | T ∈ P 0 (T ) , ∀T ∈ T n ⊂ U.
For the discrete test space V m ⊂ V, we assume that it includes the following optimal space S n := B − * P 0 (T n ) ⊂ V, i.e., V m ⊇ S n := B − * P 0 (T n ) = φ n ∈ V : B * φ n = χ n , for some χ n ∈ P 0 (T n ) .
Without any further assumptions, the following striking result show that this pair satisfies the Fortin condition. Its proof hinges on the fact that the L p duality map of any w n ∈ P 0 (T n ) is also in P 0 (T n ). Proof In view of the equivalence between the discrete inf-sup condition and the Fortin condition (see Ern & Guermond [23]), we prove this proposition by directly establishing the discrete inf-sup condition.
Let w n ∈ U n ⊆ P 0 (T n ), then Let J p (w n ) := w n 2−p p |w n | p−1 sign(w n ) denote the L p duality map of w n , and notice that it is also in P 0 (T n ). Furthermore, we have the duality-map property w n , J p (w n ) p,q = w n p J p (w n ) q . Therefore, where, in the last step, we used that B − * χ V ≤ γ −1 B χ q for all χ ∈ L q (Ω) (this is nothing but the dual counterpart of Theorem 3.A). Finally, [23,Theorem 1] implies the existence of a Fortin operator Π : V → V m with C Π = M µ /γ B .

Remark 5.4 (Petrov-Galerkin method)
If (U n , V m ) ≡ (P 0 (T n ), S n ), then dim U n = dim V m . Then, Proposition 5.3 together with (4.2b) implies that r m = 0. Thus we obtain from (4.2a) that the approximation u n satisfies the Petrov-Galerkin statement (cf. [32,Section 5]): Remark 5.5 (Cell average) If (U n , V m ) ≡ (P 0 (T n ), S n ), the approximation u n is in fact the element average of the exact solution u, i.e., To prove (5.6), note that (5.5) can be written as Let χ T be the characteristic function of the element T . Then, the test function v T = B − * χ T determines u n | T . Indeed, Remark 5.6 (Quasi-uniform meshes) In the case that U n = P 0 (T n ) where the partitions {T n } are quasi-uniform shape-regular meshes with mesh-size parameter h, the following a priori error estimate is immediate (apply Remark 4.2 with k = 0), provided that u ∈ W s,p (Ω) for 0 ≤ s ≤ 1: u − u n p h s |u| W s,p (Ω) .
It can be shown (e.g. by computing the Sobolev-Slobodeckij norm) that u ∈ W s,p (0, 1) for any 0 < s < 1/p, but not s = 1/p. Figure 4 shows the convergence of u − u n p with respect to h, for various p. The observed convergence behavior, as anticipated in Remark 5.6, is indeed close to O(h 1/p ).  Figure 4: Approximating an advection-reaction problem with a discontinuous solution using the optimal pair (P 0 (T n ), S n ): The convergence in u − u n p is close to O(h 1/p ), which is optimal for near-best approximations.
Example 5.8 (Basis for optimal test space) Let us illustrate the discrete test space S n in 1-D for the particular case where the (scalar-valued) advection β(x) is space-dependent and µ ≡ 0. Let Ω = (0, 1) and let β be a strictly decreasing and positive function such that β (x) is bounded away from zero (hence Assumption 2.1 is valid). The space V is given by Let 0 = x 1 < x 1 < · · · < x n+1 = 1 be a partition of Ω and define T n = {T j } where T j = (x j−1 , x j ). Let χ T j be the characteristic function of T j and h j = |T j |. The discrete test space S n is defined as the span of the functions v j ∈ V such that −(βv j ) = χ T j , which upon integrating over the interval [x, 1] gives : Moreover, we can combine them in order to produce the local, nodal basis functions: See Figure 5(a) for an illustration of these basis functions with h j = 0.2 for all j and β(x) = 1.001 − x.
Another interesting example is when we have two inflows, each on one side of the interval Ω = (0, 1). This is possible by means of a strictly decreasing β(x) such that β(0) > 0 and β(1) < 0, and such that β (x) is bounded away from zero.  Figure 5: Basis for the optimal test space S n that is compatible with U n = P 0 (T n ), in the case of two different space-dependent advection fields β(x) corresponding to (a) leftsided inflow and (b) two-sided inflow, respectively.
The solution u ∈ L p (Ω) of problem (3.1) may be singular at the pointx ∈ Ω for which β(x) = 0, even for smooth right hand sides. The test functions computed by solving −(βv j ) = χ T j may be discontinuous whenx matches one of the mesh points. This is illustrated in Figure 5 Example 5.9 (A practical alternative to S n ) In practise it may not be feasible to explicitly compute a basis for S n . Practical alternatives consist of, for example, continuous piecewise polynomials of sufficiently-high degree k on T n , or continuous piecewise linear polynomials on Refine (T n ), which is the submesh obtained from the original mesh T n by performing uniform refinements of all elements (see [7] for a similar alternative in a DPG setting).
In the method, we take p = 2, U n = P 0 (T n ) and V m = P 1 cont (Refine (T n )), where T n is a mesh of uniform elements of size h = 1/n, and Refine (T n ) is an -refined submesh with uniform elements of size h = h/(2 ). Figure 6 plots the convergence of the u − u n 2 versus h for = 1, 2 and 4 (error plots are actually similar for all ≥ 1). We note that = 0 is in general not sufficiently rich, as it leads to a singular matrix for h = 1/2, while the results for ≥ 1 did not show any instabilities. To anticipate the rate of convergence, note the Sobolev embedding result W s,2 (Ω) ⊂ L r (Ω) for s ≥ 1 2 − 1 r and r ≥ 2. Therefore, one expects a convergence of O(h s ) with s = 1 6 , which is indeed consistent with the to the advection-reaction problem with the DDMRes method using U n = P 0 (T n ) and V m = P 1 cont (Refine (T n )). The convergence is close to O(h 1 /6 ), which is optimal for near-best approximations.  Figure 7: Convergence of approximation u n| toward u n|∞ on a fixed mesh T n for a smooth exact solution, where u n| denotes the approximation obtained by the DDMRes method using U n = P 0 (T n ) and V m = P 1 cont (Refine (T n )). The observed convergence is O(h 2 ). numerical observation in Figure 6. The oscillations are caused by the singularity location (x = 1 /12) being closer to the left or right element edge depending on h.
To investigate for a fixed mesh with h = 1 /16 the convergence of the obtained approximations u n with respect , we consider β(x) = 2 − x, µ(x) = 0 and exact solution u(x) = 1 + 2x for x ∈ Ω. Figure 7 plots the error u n|∞ − u n| with respect to h = h/(2 ), where u n|∞ denotes the ideal approximation (V m = V). For this error we observe a rate of convergence O(h 2 ).

5.3
The pair P 0 (T n ) -P 1 conf (T n ): An optimal pair in special situations As a last application, we consider a special multi-dimensional situation such that the optimal test space S n defined in (5.4b) reduces to a convenient finite element space. We focus on a 2-D setting, and assume Ω ⊂ R 2 is polygonal and T n is a simplicial mesh (triangulation) of Ω. Let F n = {F } denote all mesh interior faces. 4 Assume that µ ≡ 0, div β ≡ 0 and that the hypothesis of Assumption 2.6 is fulfilled. Assume additionally that β is piecewise constant on some partition of Ω, and let the mesh T n be compatible with this partition, i.e., where [[·]] = (·) + − (·) − denotes the jump. Finally, assume that the mesh is flowaligned in the sense that each triangle T ∈ T n has exactly one tangential-flow face F ⊂ ∂T for which β · n T = 0 on F . Necessarily, the other two faces of T correspond to in-and out-flow on which β| T ·n T < 0 and β| T ·n T > 0, respectively. The main result for this special situation is the following characterization of S n : Proposition 5.10 (Optimal space S n : Flow-aligned case) Under the above assumptions, S n = P 1 conf (T n ) := φ n ∈ V = W q 0,+ (β; Ω) : φ n | T ∈ P 1 (T ) , ∀T ∈ T n .
Note that P 1 conf (T n ) consists of W q 0,+ (β; Ω)-conforming, piecewise-linear functions, which can be discontinuous across tangential-flow faces, but must be continuous across the other faces. Furthermore, they are zero on ∂Ω + .
Proof The proof follows upon demonstrating that B * P 1 conf (T n ) = P 0 (T n ). First note (under the above assumptions) that B * = −β · ∇ n , where ∇ n is the elementwise (or broken) gradient, i.e., (∇ n φ)| T = ∇(φ| T ) for all T ∈ T n . Since functions in P 1 conf (T n ) are element-wise linear, we thus have that B * P 1 conf (T n ) ⊂ P 0 (T n ). We next show that P 0 (T n ) ⊂ B * P 1 conf (T n ). Note that P 0 (T n ) = Span{χ T , T ∈ T n }, where χ T is the characteristic function for T . Let φ T be the unique solution in V such that B * φ T = χ T . The Ω-filling assumption (see Assumption 2.6) guarantees that β = 0 a.e. in Ω (otherwise we would have −β · ∇z ± = 0 in some element, contradicting (2.9)). Thus, for a.e. x ∈ Ω consider the polygonal path Γ(x) ⊂ Ω that starts from x and moves along the advection field β. By the Ω-filling assumption, the path Γ(x) has to end in some point on the out-flow boundary ∂Ω + (otherwise it will stays forever within Ω, contradicting the existence of a bounded function z ± ∈ W ∞ (β; Ω) whose absolute value grows linearly along Γ(x)). Hence, we can construct φ T integrating χ T over the polygonal path Γ(x) from ∂Ω + to x. By construction, φ T is a piecewise linear polynomial, which can be discontinuous only across {F ∈ F n : β · n F = 0}. Besides, φ T satisfies the homogeneous boundary condition over ∂Ω + . Hence φ T ∈ P 1 conf (T n ) and χ T ∈ B * P 1 conf (T n ). Example 5.11 (2-D numerical illustration) To illustrate the above setting with a numerical example, let Ω = (0, 1) × (0, 2) ⊂ R 2 , f • = 0, and g be nonzero on the inflow boundary ∂Ω − = {(x, 0) , x ∈ (0, 1)}. Let an initial triangulation of the domain be as in Figure 8 (top-left mesh). The advection β is such that, for the bottom, left, right and top boundary, we have that β · n is −1, 0, 0 and 1, respectively. Next, within each triangle, β is some constant vector with a positive vertical component, while satisfying the above requirements (i.e., [[β · n F ]] F = 0 on each interior face F , and each triangle has a tangential-flow, in-flow and out-flow face). 5 By Remark 5.4, the DDMRes method with spaces P 0 (T n ) and P 1 conf (T n ) can be implemented as a Petrov-Galerkin method.
We first consider the smooth inflow boundary condition g(x, 0) = sin(πx) for x ∈ (0, 1). Figure 8 (top row) shows the approximations for u n obtained on the initial triangulation and three finer meshes. The finer meshes were obtained by uniform refinements of the initial triangulation using so-called red -refinement [36, Section 2.1.2] (splitting each triangle into four similar triangles), which preserves the above flow-aligned mesh requirement. The approximations nicely illustrate the cell-average property mentioned in Remark 5.5 (the exact solution is simply found by traversing g along the characteristics). In Figure 8 (bottom) the convergence of u − u n p is shown to be optimal (rate is O(h)) for various values of p. Figure 9 shows the same results as before, but now for a discontinuous inflow boundary condition g(x, 0) = sin(πx) sign(x − 1 /3) for x ∈ (0, 1). Again the DDMRes method provides a near-best approximation; as anticipated, the observed rate of convergence is O(h 1/p ) (cf. discussion in Example 5.7).

A.1 Proof of Theorem 3.A
In this section, we give the proof of Theorem 3.A by means of the so-called Banach-Nečas-Babuška inf-sup conditions (see [22,Theorem 2.6]): Our technique is similar to the one used by Cantin [9], but note that we prove (BNB1)-(BNB2) on the adjoint bilinear form. Recall that the primal operator is a continuous bijection if and only if the adjoint operator is a continuous bijection, in which case both inf-sup constants are the same. The following proof is also analogue to the proof in Hilbert spaces given by Di Pietro & Ern [21, Section 2.1]. We start by giving some properties that we need for the Banach setting. Let J q (v) = v 2−q q |v| q−1 sign(v) ∈ L p (Ω) = U be the duality map of L q (Ω), i.e., Additionally, for any v ∈ V = W q 0,+ (β; Ω) ⊂ L q (Ω) notice the following identity: : DDMRes approximations using U n = P 0 (T n ) and V m = P 1 conf (T n ) for an incompressible advection problem with special piecewise constant β and a discontinuous inflow boundary condition g. Top row : Approximations u n on different (nested) meshes. Bottom: The convergence in u − u n p is O(h 1/p ), which is optimal for near-best approximations to the discontinuous solution u.

A.2 Proof of Proposition 5.1
We construct explicitly a Fortin operator Π : V → V m satisfying Assumption 4.1.
We note that this 1-D proof is similar to the 1-D version of the proof of [22,Lemma 4.20,p. 190]. Let −1 = x 0 < x 1 < · · · < x n = 1 be the set of nodes defining the partition T n . Over each element T j = (x j−1 , x j ) ∈ T n we define Π to be the linear interpolant Π 1 plus a quadratic bubble, i.e., x j−1 )(x − x j ). The coefficient α j multiplying the bubble Q j (x) is selected in order to fulfill the equation: Observe that Π(v) ∈ P 2 cont,0,{1} (T n ) ⊆ P k cont,0,{1} (T n ) since k ≥ 2, and for all w n ∈ U n we have: b(w n , Π(v)) = n j=1 T j w n Π(v) − w n Π(v) x j (by integration by parts) Hence, the requirement (4.4b) is satisfied. Now we recall that (·) V := (·) q . Therefore to obtain the requirement (4.4a) (i.e. the boundedness of the operator Π), we note that on each element: Thus, on each element (and therefore globally) we have: where the constant C q = 6/(q + 1) 1 q is mesh-independent.