SUPG-stabilized Virtual Elements for diffusion-convection problems: a robustness analysis

The objective of this contribution is to develop a convergence analysis for SUPG-stabilized Virtual Element Methods in diffusion-convection problems that is robust also in the convection dominated regime. For the original method introduced in [Benedetto et al, CMAME 2016] we are able to show an"almost uniform"error bound (in the sense that the unique term that depends in an unfavorable way on the parameters is damped by a higher order mesh-size multiplicative factor). We also introduce a novel discretization of the convection term that allows us to develop error estimates that are fully robust in the convection dominated cases. We finally present some numerical result.


Introduction
The Virtual Element Method (VEM) was introduced in [6,7] as a generalization of the Finite Element Method (FEM) to general polygonal and polyhedral meshes. Since its introduction, the VEM enjoyed a wide success in the Numerical Analysis and Engineering communities, both due to the encouraging results and the natural construction.
The possibility of using general polytopal meshes makes VEM suitable for diffusion problems, for instance by making much easier to adapt to complex geometry of the data (such as in basin and reservoir simulations) and to irregularities of the solution. The VEM literature on the diffusion-reaction-convection problem is indeed very wide, covering primal and mixed methods, conforming and non-conforming schemes, ranging from foundation/theoretical contributions to more applicative articles; a very short representative list being [9,8,20,5,24,10,16,14,30,31,12,26]. Some examples of other numerical methods for the diffusion-reactionconvection problem that can handle polytopal meshes are [27,28,4,22,3,21]. On the other hand, the majority of the VEM contributions assume a dominant diffusion and do not address the significant case of convection dominated problems. Indeed, as it happens for standard FEM, unless some ad-hoc modification is introduced, also the VEM is expected to suffer in convection dominated regimes, leading to very large errors unless the mesh is extremely fine. To the best of the authors' knowledge, only in the papers [15,17] such issue is addressed; in these articles a SUPG-stabilized Virtual Element scheme for conforming and non-conforming VEM is proposed, and analyzed both theoretically and numerically. However, the stability and convergence analysis in [15,17] is not uniform in the diffusion/convection parameters, and therefore it cannot be used to theoretically justify the method behaviour in the convection dominated regime. Moreover, a sufficiently small mesh size h is required to carry out the analysis. The main difficulty in deriving uniform error estimates for SUPG-stabilized VEM is handling a variable convection coefficient in the presence of projection operators (which are needed in the VEM construction), that partially disrupt the structure of the convection term.
The aim of the present paper is to address, in the conforming case, this challenging theoretical aspect, thus deriving convergence estimates for a slight modification of the SUPG VEM scheme of [15] that are robust in the involved parameters and do not require a sufficiently small h condition. We think that, in addition to filling an important theoretical gap, having this deeper understanding is fundamental in order to develop SUPG stabilizations in more complex settings, such as fluid-dynamics problems. For instance, deriving the aforementioned proofs inspired us to propose also a novel (alternative) approach for the discretization of the convective term, in addition to the original one. For the (slightly modified) discrete convection form introduced in [15], we are able to show an error estimate that is "almost uniform" in the involved parameters, in the sense that the unique term that depends in an unfavourable way on the parameters is damped by a higher order multiplicative factor in h. For the novel form here proposed, we are able to show full robustness in the parameters. Finally, for the sake of completeness we also present a few numerical results, the main objective being to make a practical comparison among some different discretization options described in the previous section.
The present paper is organized as follows. In Section 2 we present the continuous problem and in Section 3 we introduce some preliminaries and notation. Afterwards, in Section 4 we review the SUPG-stabilized Virtual Element Method under analysis, also introducing the novel convective term option. The main contribution of this article is Section 5, where we develop the aforementioned convergence analysis. The numerical tests are shown in Section 6.
Throughout the paper, we will follow the usual notation for Sobolev spaces and norms [1]. Hence, for an open bounded domain ω, the norms in the spaces W s p (ω) and L p (ω) are denoted by · W s p (ω) and · L p (ω) respectively. Norm and seminorm in H s (ω) are denoted respectively by · s,ω and |·| s,ω , while (·, ·) ω and · ω denote the L 2 -inner product and the L 2 -norm (the subscript ω may be omitted when ω is the whole computational domain Ω).

Continuous Problem
Let Ω ⊂ R 2 be the computational domain and let ε > 0 represent the diffusive coefficient (assumed to be constant), while β ∈ [L ∞ (Ω)] 2 with divβ = 0, is the transport advective field, and f ∈ L 2 (Ω) is the volume source term. Then, our linear steady advection-diffusion model problem reads where V = H 1 0 (Ω) and the bilinear forms a(·, ·) : By a direct computation, being divβ = 0, it is easy to see that the bilinear form b(·, ·) is skew Therefore, the bilinear form b(·, ·) is equal to its skew-symmetric part, defined as: However, at the discrete level b(·, ·) and b skew (·, ·) will lead to different bilinear forms, in general. It is well known that discretizing problem (1) leads to instabilities when the convective term β [L ∞ (Ω)] 2 is dominant with respect to the diffusive term ε (see for instance [33]). In such situations a stabilized form of the problem is required in order to prevent spurious oscillations that can completely spoil the numerical solution. In the following sections we propose a virtual elements version of the classical Streamline Upwind Petrov Galerkin (SUPG) approach [32,29]. From now on, we assume that the material parameters are scaled so that it holds: We finally remark that the proposed approach can be trivially extended to more general situations such as reaction-convection-diffusion problems, non-constant diffusive coefficients and different boundary conditions. We here prefer to focus on the fundamental difficulties of the problem rather than deploying a "smokescreen" of additional minor technicalities. Moreover, also the analysis of the three dimensional case could be developed with very similar arguments.

SUPG stabilizing form
From now on, we will denote with E a general polygon, e will denote a general edge of E, moreover |E| and h E will denote the area and the diameter of E respectively, whereas n E will denote the unit outward normal vector to ∂E.
We suppose that { Ω h } h fulfils the following assumption: (A1) Mesh assumption. There exists a positive constant such that for any • any edge e of E has length ≥ h E .
We remark that the hypotheses above, though not too restrictive in many practical cases, could possibly be further relaxed, combining the present analysis with the studies in [11,19,25,13].
We now briefly review the construction of the SUPG stabilization [32,29] for the advectiondominated problem (1). First of all, we decompose the bilinear forms a(·, ·) and b skew (·, ·) into local contributions, by defining The global approximated bilinear form A supg (·, ·) and the global right-hand side are defined by simply summing the local contributions: Since the exact solution u of equation (1) is well defined for all v ∈ V and u solves the stabilized problem The aim of the following sections is to derive a VEM discretization of the stabilized problem (11). In the following the symbol will denote a bound up to a generic positive constant, independent of the mesh size h, of the SUPG parameter τ E , of the diffusive coefficient ε and of the transport advective field β, but which may depend on Ω, on the "polynomial" order of the method k and on the regularity constant appearing in the mesh assumption (A1).

Projections and polynomial approximation properties
In the present subsection we introduce some basic tools and notations useful in the construction and the theoretical analysis of Virtual Element Methods.
Lemma 3.2 (Inverse estimate). Let, for any E ∈ Ω h , γ E denote the smallest positive constant such that for any p n ∈ [P n (E)] 2 , it holds Then, under assumption (A1), there exists γ ∈ R + such that

Virtual Element spaces
Let k ≥ 1 be the "polynomial" order of the method. For any E ∈ Ω h we consider the local "enhanced" virtual element space [2] given by (15) We here summarize the main properties of the space V h (E) (we refer to [2] for a deeper analysis).
(P3) Polynomial projections: the DoFs D V allow us to compute the following linear operators: The global virtual element space is obtained by gluing such local spaces, i.e.
with the associated set of degrees of freedom. We finally recall from [23,19] the optimal approximation property for the space V h (Ω h ).

Virtual Element Forms
The next step in the construction of our method is to define a discrete versions of the stabilized SUPG form A supg (·, ·) in (5). It is clear that for an arbitrary pair (u is not computable since u h and v h are not known in closed form. Therefore, following the usual procedure in the VEM setting, we need to construct a computable discrete bilinear form. In the following, in accordance with definition (5), we define a discrete counterpart of each brick composing A E supg .
Here, the stabilizing bilinear form for two positive uniform constants α * and α * . The condition above essentially requires the For instance, the standard choices for the stabilization are the dofi-dofi stabilization [6] and the D-recipe stabilization introduced in [10].
Concerning the approximation of the convective term b E (·, ·), we here propose two possible choices: recalling property (P3), let us define for all While the form (19) follows a more standard "approximation by projection" VEM approach (see for instance [15]), the novel form (20) is amenable to the development of an improved theoretical result. In the following b E h (·, ·) : V h (E) × V h (E) → R will denote indifferently one of the aforementioned forms and, in accordance with (4) Exploiting again property (P3), the stabilized forms B E (·, ·) in (6) and L E (·, ·) in (7) are discretized as follows where β E := β [L ∞ (E)] 2 and the parameter τ E > 0 has to be chosen. In accordance with (5), the VEM stabilized form The corresponding computable VEM version of the SUPG right-hand side in (8) reads as (26) and its global counterpart is

Virtual Element SUPG problem
Referring to the discrete space (16), the discrete bilinear form (25) and the approximated right-hand side (27), the virtual element SUPG approximation of the advection-dominated diffusion equation (1) is

Theoretical analysis
In this section we analyze the stabilization method defined in (28). In particular, we assess the stability property of problem (28) and we provide the convergence error estimate for the discrete solution obtained with both discrete convective forms defined in (19) and (20). All estimates clearly display the dependence on the mesh size h, on the parameter τ E and the problem data ε and β.

Stability
Let us start with the stability analysis for the proposed VEM SUPG method. First of all we define the VEM SUPG norm

Proposition 5.1 (Coercivity). Under the assumption (A1) if the parameters τ E satisfy
where γ E is the constant appearing in the inverse estimate of Lemma 3.2, the bilinear form Proof. We simply consider all the terms in the sum (24). For the first three terms by definitions (17), (21) and (22) and stability estimate (18) we get whereas for the last term we infer Collecting the previous bound, (32) and (33) we obtain Remark 5.1. Notice that the norm · supg,E is slightly different from the usual norm introduced in standard SUPG theory [32,29] However we observe that the "classical norm" · 2 supg,E is controlled by the "VEM norm" · supg,E . Indeed, recalling (34), for any

Error estimates
The aim of the present section is to derive the rate of convergence for the proposed SUPG virtual element scheme (28) in terms of the mesh size h, the SUPG parameter τ E , the diffusive coefficient ε and transport advective field β. The hidden constants may depend on Ω, on k, on the regularity constant appearing in the mesh assumption (A1) and on the stability constants α * and α * (cf. (18)).
Let u ∈ V and u h ∈ V h (Ω h ) be the solutions of problem (11) and problem (28), respectively, and let us define the following error functions is the interpolant function of u defined in Lemma 4.1, and Π ∇ k u ∈ P k (Ω h ) is the piecewise polynomial defined in (14). We introduce the analysis with the following abstract error estimation.
The next step in the analysis consists in estimating all the terms in the bound (35). We make the following assumption: (A2) Data assumption. The solution u, the advective field β and the load f in (11) satisfy: for some 0 < s ≤ k.
Note that in the following lemmas it is not restrictive to assume β E > 0 since β E = 0 implies β| E = 0 and thus the corresponding terms vanish. can be bounded as follows (for 0 < s ≤ k) Proof. Applying the definition of the norm · supg,E , of the L 2 -orthogonal projection Π 0,E k−1 (cf. (12)), of the H 1 -orthogonal projection Π ∇,E k (cf. (13)), and the interpolation estimate of Lemma 4.1, we easily obtain The thesis now follows by summing the local contributions.
. Under the assumptions (A1) and (A2), the term η E F can be bounded as follows (for 0 < s ≤ k) where for any E ∈ Ω h λ E := min 1 Proof. Applying the definitions (8), (26) and the definition of L 2 -orthogonal projection we obtain Using a scaled Poincaré inequality we infer Recalling the definition of the norm · supg,E and the stability of Π ∇,E k with respect to the H 1 -seminorm, from Lemma 3.1 we get Regarding the second term η E F ,2 , from (34) and Lemma 3.1 we obtain Now the thesis follows from (36), (37) and (38). 2 represents a locally scaled regularity term for β. Roughly speaking, it is related to the local variations of β and not to its amplitude.

Lemma 5.3 (Estimate of η E a ). Under the assumptions (A1) and (A2), the term η E a can be bounded as follows (for
Proof. The consistency and the continuity of the form a E h (·, ·), Lemma 3.1 and Lemma 4.1 easily imply (A2), the term η E B can be bounded as follows (for 0 < s ≤ k)

Lemma 5.4 (Estimate of η E B ). Under the assumptions (A1) and
Proof. Using the definition of L 2 -projection, simple computations yield We analyse separately each term in the sum. The term η E B,1 is bounded using (34) and the continuity of Π 0,E k−1 with respect to the L 2 -norm, Lemma 3.1 and Lemma 4.1: For the second term η E B,2 using again (34) and Lemma 3.1 we infer Finally for the last term in (39), employing (18), the stability of the H 1 -seminorm projection with respect to the H 1 -seminorm, Lemma 3.1 and Lemma 4.1 we get The thesis now follows by collecting (40), (41) and (42) in (39). (A2), the term η E L can be bounded as follows (for 0 < s ≤ k)

Lemma 5.5 (Estimate of η E L ). Under the assumptions (A1) and
Proof. By definition of L 2 -orthogonal projection we infer The term η E L,1 , employing Lemma 3.2, Lemma 3.1 and Lemma 4.1 is estimated as follows The second term in (43), recalling (34), can be easily bounded as follows Proof. By the definition of the skew symmetric forms (4) and (21) we need to estimate the terms We now analyse each term η E b,i for i = 1, . . . , 6 in the sum above. • η E b,1 : using an integration by parts, bound (34) and the definition of · supg,E we infer (49) • η E b,2 : a scaled Poincaré inequality and the definition of L 2 -projection imply • η E b,3 : from the definition of L 2 -projection, the Poincaré inequality and Lemma 3.1, we infer • η E b,4 : using similar computations of the previous item we obtain • η E b,5 : exploiting the property of L 2 -projection and bound (34) we get • η E b,6 : using similar computations of the previous item we have (54) The thesis now follows gathering (49)-(54) in (48).
be the bilinear form in (20). Then under assumptions (A1) and (A2), the term η E b,∂ can be bounded as follows where σ E is defined in (47).
Proof. Recalling definition (20) we need to estimate the terms By integration by parts we have We now analyse each term η E b,i for i = 1, . . . , 4 in the sum above. • η E b,1 : using the same computations in (50) we infer • η E b,2 : exploiting the computation in (53) we obtain • η E b,3 + η E b,4 : we use a scaled trace inequality [18] making use of the scaled norm |||v|||  (19), we are able to obtain a "damped" dependence on ε: in estimate (46) the term ε −1/2 |β| [W 1 ∞ (E)] 2 h s+2 E |u| s+1,E blows up as ε → 0, but at the same time it is of higher order with respect to h E . Instead, for the new form (20) we are able to obtain full independence from ε.
We are now ready to prove the convergence results for the proposed VEM SUPG scheme. The error estimates in Lemmas 5.1-5.7 are explicit in the parameters of interest: the mesh size h, the diffusive coefficient ε, the advective field β and the SUPG parameter τ E . In order to simplify the final estimate and to make clearer the implications of the convergence results, in the following propositions we include the Sobolev regularity terms for u, f and the normalized Proof. The proof is a direct consequence of Proposition 5.2, Lemmas 5.1, 5.3, 5.4, 5.5, and 5.6 we estimated the last two terms in (46) as follows. The penultimate term is bounded using the Poincaré inequality on the domain Ω For the last term, noticing that e I , e h ∈ V , it holds that  (20). Then it holds that Proof. The proof follows from Proposition 5.2, Lemmas 5.1, 5.3, 5.4, 5.5, 5.7 and equation (59), also recalling that We conclude that in the diffusion dominated regime both schemes yield the optimal rate of convergence. In the convection dominated regime only the scheme derived from the bilinear forms b E ∂,h (·, ·) has the optimal rate of convergence. For the scheme derived from b E o,h (·, ·) the error is polluted by ε −1 . Nevertheless, we stress that such a factor appears in front of the "higher" order term h 3 , therefore the influence of the diffusion coefficient is strongly reduced.

Numerical experiments
In this section we numerically validate the proposed methods by means of the following model problem.
We choose the boundary conditions and the source term (which turns out to depend on ε) in such a way that the analytical solution is always the function u(x, y) := sin(π x) sin(π y) .
Guided by the definition of the || · || supg norm (cf. (29) and (30)), by the error estimates of Propositions 5. 3-5.4, and noticing that the discrete solution u h ∈ V h (Ω h ) is not explicitly pointwise available, the following error quantities will be considered.
• H 1 −seminorm error • convective norm error As far as the mesh types are concerned, we take the following: • quad: a mesh composed by structured quadrilaterals; • tria: a Delaunay triangulation of the unit square; • voro: a centroidal Voronoi tessellation of the unit square where the cells' shape is optimized via a Lloyd algorithm • rand: a voronoi tessellation of the unit square where the cell shapes are not optimized.
In Figure 1 we show an example of such meshes, and we also remark that the former two types can be used in connection with a standard finite element procedure, contrary to the latter two. Moreover, since we are interested in a robustness analysis with respect to the diffusion parameter, for each mesh sequence we take ε = 10 −3 and 10 −6 .
Effect of the SUPG stabilization. Before assessing the convergence properties of the proposed methods, we check the effect of inserting the SUPG term in the variational formulation of the problem. Here we focus on the first choice for the convective bilinear form, without the algebraical skew-symmetrization; namely we here use form (19). In Figures 2 and 3 we show the convergence graphs for different approximation degrees k and for each type of the meshes. The cases with and without the SUPG term are labelled as SUPG and NONE, respectively. As expected, dropping the SUPG term clearly deteriorates the quality of the discrete solutions for all approximstion degree k, the degeneration being heavier for smaller ε. We do not report the results for "large" ε, but we have experienced that, of course, both the SUPG and the NONE approaches give similar outcomes, attaining the correct convergence rate.
Incidentally, we remark that Figures 2 and 3 also display the robustness of the present SUPG virtual element approach with respect to element shape and distortion. In fact, given an approximation degree k, the convergence histories are rather similar for all the mesh types, even though the rand meshes contain some odd shaped elements. Different discrete convective bilinear forms. We now present the numerical results for different variants of the convective term approximation. More precisely, we consider four types of (local) discrete convective terms: b E o,h (·, ·) and b E ∂,h (·, ·), and their skew-symmetric counterparts, see (19), (20) and (21). We refer to such approximations as orig and boun, respectively, while we label their skew-symmetric counterparts as origSkew and bounSkew (the ones covered by the theoretical analysis developed in Section 5).
To avoid instabilities in the convection dominated regime, in what follows we always take advantage of the SUPG stabilization.
In Figures 4 and 5 we collect the results with the diffusion coefficients ε = 10 −3 and ε = 10 −6 , respectively. We notice that the convergence rate of both error norms e C and e H 1 are the expected ones, and they are robust in the parameter ε. When we consider the H 1 -seminorm error e H 1 , we observe that the origSkew behaves slightly worst than the other three discrete forms, especially for the approximation degrees k = 2, 3, and independently of the values of ε.
We conclude the section by noticing that b E ∂,h (·, ·) and b E o,h (·, ·) coincide for a costant convection term β (the proof can be easily performed by a direct computation). To have a numerical evidence about this fact, we have considered a problem with constant vector field β and compared the stiffness matrices provided by b E o,h (·, ·) and b E ∂,h (·, ·), respectively. For every mesh and every approximation degree, we have found that they always differ up to machine precision. Here we show the data (norms of the difference between the stifness matrices) only for the finest voro mesh and for k = 1, 2, see Table 1.