SUPG-stabilized stabilization-free VEM: a numerical investigation

1.   Introduction

The development of numerical methods for the solution of partial differential equations exploiting general polygonal or polyhedral meshes has been a topic of great interest in later years. Among the many families of schemes developed in this context ^[1,2,3,4,5], this paper considers Virtual Element Methods (VEM).

Since the seminal papers ^[6,7,8,9], VEM have been applied to many problems in which polytopal meshes can be exploited in order to better handle the geometrical complexities of the computational domain. A brief, not exhaustive, list of papers is ^{[10,11,12,13,14,15,16,17,18,19]}. Virtual element schemes are based on the definition of locally computable polynomial projections that are involved in the discrete bilinear forms. These forms consist of the sum of two terms: a singular one that is consistent on polynomials and a stabilizing one that ensures coercivity. In the literature, the arbitrary nature of the stabilization term remains an issue to be investigated. Indeed, it has been shown that it can cause problems in many theoretical and numerical contexts. For instance, the isotropic nature of the stabilization can become an issue when devising Streamline Upwind Petrov-Galerkin (SUPG) stabilizations ^[20,21].

In the context of the theoretical development of a VEM scheme that does not require a stabilization term, Stabilization-Free Virtual Elements Methods (SFVEM) have been recently introduced in ^[22,23], for the lowest order primal and mixed discretizations of the Poisson equation. The key idea of the proposed method is to define self-stabilized bilinear forms, exploiting only higher order polynomial projections. The theoretical study of this method is an ongoing investigation, however some tests on highly anisotropic problems show that this new approach is able to overcome some issues of standard VEM related to the isotropic nature of the stabilization ^[24,25]. Moreover, this approach has been successfully applied to linear and non-linear elasticity problems in ^{[26,27,28,29]}. The aim of this paper is to numerically investigate the SFVEM in the context of the SUPG formulation of advection-diffusion problems in the advection-dominated regime.

The outline of the paper is as follows: In Section 2 we present the advection-diffusion model problem. In Section 3 we define the numerical scheme. In Section 4 we discuss the well-posedness of the discrete problem. Section 5 is devoted to a priori error estimates. Finally, in Section 6 we present some numerical results that assess the behaviour of the proposed method, also comparing it to standard SUPG-stabilized VEM schemes.

In the following, $\left({\cdot}, {\cdot}\right)_{\omega}$ denotes the $\mathrm{L}^{2}{\!\left({\omega}\right)}$-scalar product, $\left\lVert{{\cdot}}\right\rVert_{\omega}$ denotes the corresponding norm, and $\left\lVert{{\cdot}}\right\rVert_{m, \omega}$ and $\left\lvert{{\cdot}}\right\rvert_{m, \omega}$ denote the $\mathrm{H}^{{m}}_{}{\!({\omega})}$ norm and semi-norm.

2.   Model problem

Let $\Omega \subset \mathbb{R}^2$ be a bounded open set. We consider the following advection-diffusion model problem: find $u$ such that

where $\varepsilon_{} > 0$ is a positive real number and we assume that $\beta \in \left[\mathrm{L}^{\infty} {\!\left({\Omega}\right)}\right]^2$, $\nabla\cdot\beta = 0$, and $f\in\mathrm{L}^{2} {\!\left({\Omega}\right)}$. Moreover, we define the bilinear form ${{B}}\colon \mathrm{H}^{{1}}_{0} {\!({\Omega})}\times \mathrm{H}^{{1}}_{0} {\!({\Omega})}\to \mathbb{R}$ and the operator $F^{}\colon \mathrm{L}^{2} {\!\left({\Omega}\right)}\to \mathbb{R}$ such that

Then, the variational formulation of (2.1) reads as follows: find $u\in \mathrm{H}^{{1}}_{0} {\!({\Omega})}$ such that

It is a standard result that the above problem (2.2) is well-posed under the above regularity assumptions on the data. In the following, we consider homogeneous Dirichlet boundary conditions, but more general boundary conditions can be considered and will be considered in the numerical tests. Finally, for any open set $\omega \subseteq \Omega$, we define:

3.   Problem discretization

This section is devoted to the discretization of (2.2) using an enlarged enhancement Virtual Element space. The discretization defined here was introduced in ^[24].

3.1.   Discrete space

Let $\mathcal{M}_{h }$ be a polygonal mesh of $\Omega$. Let $h_E$ denote the diameter of $E\in\mathcal{M}_{h }$, and let $h = \max_{E\in\mathcal{M}_{h }} h_E$ be the mesh parameter. We make the following mesh assumptions ^[6,9]: There exists a constant $\kappa > 0$ independent of $h$ such that, for any polygon $E\in\mathcal{M}_{h }$, if ${\mathcal{E}_{{h}}^{E}}$ denotes the set of edges of $E$,

$(1)$ $E$ is star-shaped with respect to a ball of radius $\rho \geq \kappa h_E$;

$(2)$ $\forall e \in{\mathcal{E}_{{h}}^{E}}$, $\left\lvert{{e}}\right\rvert \geq \kappa h_E$, where $\left\lvert{{e}}\right\rvert$ denotes the length of $e$.

For each $E\in\mathcal{M}_{h }$, let $\mathbb{P}_{{n}}\!\left({E}\right)$ be the space of polynomials of degree at most $n$, for any $n\in\mathbb{N}$. As a basis of $\mathbb{P}_{{n}}\!\left({E}\right)$, we choose the following set of scaled monomials:

where $(x_E, y_E)$ is the center with respect to which $E$ is star-shaped. Moreover, let us define

Let $\Pi^{\nabla, E}_{{n}}\colon \mathrm{H}^{{1}}_{} {\!({E})}\to \mathbb{P}_{{n}}\!\left({E}\right)$ be the projection operator such that, $\forall v \in \mathrm{H}^{{1}}_{} {\!({E})}$,

Then, given $\ell_E\in\mathbb{N}$ we introduce the enlarged enhancement VEM discrete space:

where $\gamma$ denotes the trace operator and

Remark 3.1. The space $V_{{k, 0}}^{E}$ is the standard VEM space used in ^[9,20].

Notice that, $\forall \ell_E$, a possible set of degrees of freedom of $V_{{k, \ell_E}}^{E}$ is the one used for standard VEM spaces (see e.g., ^[9,20]), that is, for any $v_h\in V_{{k, \ell_E}}^{E}$,

$(1)$ the values of $v_{h}$ at the vertices of $E$;

$(2)$ for each $e\in{\mathcal{E}_{{h}}^{E}}$, the values of $v_{h}$ at $k-1$ points internal to $e$;

$(3)$ the scaled moments $\frac{1}{\left\lvert{{E}}\right\rvert}\left({v_{h}}, {m_{\boldsymbol{\alpha}}}\right)_{E}$, $\forall m_{\boldsymbol{\alpha}}\in\mathcal{M}_{{k-2}}\!\left({E}\right)$.

3.2.   Discrete problem

This section is devoted to define the SUPG-stabilized discrete formulation of (2.2), following ^[20].

Let $E\in \mathcal{M}_{h }$ be given. If $k > 1$, let $\tilde{C}_{k}$ be the largest constant, independent of $h_E$, such that

where $\tilde{C}_{k}$ depends on the constant $\kappa$ appearing in the mesh assumptions, and on the degree $k$ (see e.g., ^[30]). The Péclet number associated to $E$ is defined as

Following the standard SUPG approach ^[31], for any $E\in\mathcal{M}_{h }$ let the parameter $\tau_E$ be given, such that

To define the SUPG formulation of the problem, we introduce the space

and $\forall E \in \mathcal{M}_{h }$, the bilinear form ${{{\mathcal{B}}}^{{E}}_{{}}} \colon \mathrm{H}^{{1, \Delta}}_{\mathrm{loc}} {\!({{\mathcal{M}_{h }}})}\times \mathrm{H}^{{1}}_{0} {\!({\Omega})}\to \mathbb{R}$ such that, $\forall w \in \mathrm{H}^{{1, \Delta}}_{\mathrm{loc}} {\!({{\mathcal{M}_{h }}})}$, $\forall v \in \mathrm{H}^{{1}}_{0} {\!({\Omega})}$,

where

Moreover, we introduce the operator $\mathcal{F}_{}^{E}\colon \mathrm{H}^{{1}}_{} {\!({E})}\to \mathbb{R}$ such that

The above operators are not computable from the degrees of freedom for functions in $V_{{k, \ell_E}}^{E}$. For this reason, we define discrete bilinear forms involving polynomial projections. For a given $n\in\mathbb{N}$ let $\Pi^{0, E}_{{n}}\colon \left[\mathrm{L}^{2} {\!\left({{E}}\right)}\right]^2\to \left[\mathbb{P}_{{n}}\!\left({E}\right)\right]^2$ denotes the element-wise $\mathrm{L}^{2} {\!\left({E}\right)}$ projection onto $\mathbb{P}_{{n}}\!\left({E}\right)$ and let

Remark 3.2. Notice that, contrarily to standard VEM SUPG formulations (cf. ^[20,21]), here we do not require an additional non-polynomial stabilization term in the bilinear form ${{{a}}^{{E}}_{{h}}}$. The choice of the parameter $\ell_E$ is done in order to guarantee coercivity of ${{{a}}^{{E}}_{{h}}}$, as detailed in the next section.

Finally, to state our discrete problem, we define the bilinear form ${{{\mathcal{B}}}^{{E}}_{{h}}} \colon \mathrm{H}^{{1, \Delta}}_{\mathrm{loc}} {\!({{\mathcal{M}_{h }}})}\times \mathrm{H}^{{1}}_{0} {\!({\Omega})}\to \mathbb{R}$ such that, $\forall w \in \mathrm{H}^{{1, \Delta}}_{\mathrm{loc}} {\!({{\mathcal{M}_{h }}})}$, $\forall v \in \mathrm{H}^{{1}}_{0} {\!({\Omega})}$,

and the operator $\mathcal{F}_{h}^{E}\colon \mathrm{H}^{{1}}_{} {\!({E})}\to \mathbb{R}$ such that

Then, given

our discretization of (2.2) reads: find $u_h \in V_{{k, \ell}}^{}$ such that

4.   Well-posedness

The aim of this section is to discuss the key points needed for the well-posedness of (3.5). It is currently an open problem to identify a robust criterium to choose $\ell_E$ for any kind of polygon, in such a way that ${{{a}}^{{E}}_{{h}}}$ is coercive. Thus, we assume that there exists at least a good choice of $\ell_E$ for each type of polygon and in Section 6 we perform a numerical investigation of it.

Assumption 4.1. We assume that, $\forall k \geq 1$ $\exists \ell_E$ such that $\exists \alpha > 0$ independent of $h_E$ satisfying

From now on, we set $\ell_E$ as the smallest integer satisfying Assumption 4.1. With this choice, we define the following VEM-SUPG norm:

Then, we can prove the following well-posedness result.

Theorem 4.1. Under Assumption 4.1, we have, $\forall E\in \mathcal{M}_{h }$ and for $h$ sufficiently small,

Proof. Let $v_h\in V_{{k, \ell}}^{}$ be given. First, exploiting the definition of $\tau_E$ in (3.4), the inverse inequality (3.3) and Young inequalities, we get

Thus, it follows

Moreover, since ${{{b}}^{{}}_{{}}} {\left({{v_h}}, {{v_h}}\right)} = 0$, we get

where $C_{\varepsilon_{}\beta}$ depends in particular on the problem data and the equivalence constant in (4.1). Collecting the latest estimate and (4.2), we get the thesis:

□

Notice that the above result is analogous to the one obtained in ^[20] in the case of standard VEM. Other choices are possible to discretize the transport term (see ^[32]), and they would lead to a similar well-posedness results.

5.   Error analysis

In this section we address the a priori error analysis of the proposed scheme, under Assumption 4.1. The analysis follows the techniques already used to prove a priori estimates for standard VEM schemes ^[9,20,32].

We provide some details on the interpolation properties of the considered space, since it is a slight modification of the classical VEM spaces, due to the enlarged enhancement property of $V_{{k, \ell_E}}^{E}$. The proof follows the one in ^[33] for the interpolation on standard VEM spaces. First, we prove an auxiliary result, introducing an inverse inequality in $V_{{k, \ell_E}}^{E}$.

Lemma 5.1. Let $E\in\mathcal{M}_{h }$ and let $w\in\mathrm{H}^{{1}}_{} {\!({E})}$ such that $\Delta w\in\mathbb{P}_{{k+\ell_E}}\!\left({E}\right)$, for some chosen $\ell_E \geq 0$. Then, there exists a constant $C_I$ such that

The constant $C_I$ depends on $k$, $\ell_E$ and on the mesh regularity parameter $\kappa$.

Proof. First, let $p\in\mathbb{P}_{{k+\ell_E}}\!\left({E}\right)$. Then, let $\psi_E \in\mathrm{H}^{{1}}_{0} {\!({E})}$ be a bubble function defined on a regular sub-triangulation of $E$ (see ^[33]), such that $1\geq\psi_E\geq 0$. Then, since $\psi_E p\in \mathrm{H}^{{1}}_{0} {\!({E})}$,

where we use the fact that, since $\psi_E$ is non-negative and bounded, $\sqrt{\left({\psi_E}, {p^2}\right)_{E}}$ is a norm on $\mathbb{P}_{{k+\ell_E}}\!\left({E}\right)$, and thus equivalent to $\left\lVert{{p}}\right\rVert_{E}$. Similarly, since $\psi_E\in \mathrm{H}^{{1}}_{0} {\!({E})}$, $\left\lVert{{\nabla(\psi_Ep)}}\right\rVert_{E}$ is a norm on $\mathbb{P}_{{k+\ell_E}}\!\left({E}\right)$, and standard scaling arguments provide the weight $h_E^{-1}$. It follows that $C_B$ depends on $k$, $\ell_E$ and on the shape regularity of $E$. Taking $p = \Delta w$ in the above result and applying Green's theorem and a Cauchy-Schwarz inequality, we conclude, defining $C_I = C_B^{-1}$,

□

Theorem 5.1. Let $u \in \mathrm{H}^{{s}}_{} {\!({\Omega})}$, $1\leq s \leq k+1$. Let $\ell_E \in \mathbb{N}$ be given $\forall E \in \mathcal{M}_{h }$. There $\exists u_I\in V_{{k, \ell}}^{}$ such that

Proof. Following ^[33], let $\mathcal{T}_{h }$ be the sub-triangulation of $\mathcal{M}_{h }$ obtained as the union of local sub-triangulations of each polygon $E\in\mathcal{M}_{h }$, linking each vertex to the center of the ball with respect to which $E$ is star-shaped. $\mathcal{T}_{h }$ inherits the shape-regularity of $\mathcal{M}_{h }$. Let $u_C \in \mathbb{P}_{{k}}\!\left({\mathcal{T}_{h }}\right)$ be the piecewise polynomial Clément interpolant of $u$ over $\mathcal{T}_{h }$. It holds true that (see [34, Theorem 1])

where $C_{Cl, k}$ depends on the shape-regularity of $\mathcal{M}_{h }$ and the order $k$. Let $w_I\in \mathrm{H}^{{1}}_{} {\!({\Omega})}$ be the function that solves, $\forall E \in \mathcal{M}_{h }$,

By the definition of $w_I$ we have that, $\forall E \in \mathcal{M}_{h }$, $\Pi^{0, E}_{{k}}u_C - w_I$ solves the following problem:

It follows that

which implies, exploiting the continuity of the operator $\Pi^{0, E}_{{k}}$ and (5.3),

where $C_{\Pi, k}$ depends on the shape-regularity of $E$ and the order $k$. We define $u_I \in V_{{k, \ell}}^{}$ such that $\forall E \in \mathcal{M}_{h }$

Notice that applying Green's theorem we get from (5.6) and (5.7) that $\Pi^{\nabla, E}_{{k}}u_I = \Pi^{\nabla, E}_{{k}}w_I$. Thus (5.8) implies

We now prove (5.2) for $u_I$. Concerning $\mathrm{H}^{{1}}_{} {\!({\Omega})}$-seminorm of $u-u_I$, recalling (5.3) we get

Moreover, from the definition of $w_I$ (5.4) and the definition of $u_I$ (5.6) we get $\gamma_{e} {}{\!\left({u_C}\right)} = \gamma_{e} {}{\!\left({w_I}\right)} = \gamma_{e} {}{\!\left({u_I}\right)}$ for each edge $e$ of $\mathcal{M}_{h }$. Thus,

We can thus estimate the $\mathrm{L}^{2} {\!\left({\Omega}\right)}$-norm of $u-u_I$ as follows:

where $C_{p, E}$ depends on the shape-regularity of $E$. We now focus on estimating $\left\lVert{{\nabla(u_C-u_I)}}\right\rVert_{E}$ $\forall E \in \mathcal{M}_{h }$. Applying (5.5) we get

Finally, the last term can be bounded applying Green's theorem (recall that $w_I-u_I\in\mathrm{H}^{{1}}_{0} {\!({E})}$), (5.7), (5.9), the fact that $\Pi^{\nabla, E}_{{k}}u_I = \Pi^{\nabla, E}_{{k}}w_I$, a Cauchy-Schwarz inequality, the approximation properties of polynomial projections, (5.5) and the inverse inequality (5.1) we get

Then, collecting (5.10), (5.11) and the last estimate, we obtain (5.2). □

The interpolation estimate provided by Theorem 5.1 along with approximation results analogous to the ones obtained in ^[20,32] are used to prove the following a priori error estimates, whose proof is omitted since it is analogous to the one in the cited references.

Theorem 5.2. Assume $u\in\mathrm{H}^{{s+1}}_{} {\!({\Omega})}$ and $f\in\mathrm{H}^{{s-1}}_{} {\!({\Omega})}$. Then, under the current regularity assumptions and if $\ell_E$ is chosen $\forall E\in\mathcal{M}_{h }$ in such a way that (4.1) holds,

where $C$ is independent of $h$ and on the problem coefficients and $C_{f, \varepsilon_{}\beta}$ depends on local variations of the problem coefficients.

6.   Numerical results

In this section we present three numerical experiments. In the first and second tests we confirm the convergence rates predicted by the a priori error analysis of Section 5 and we compare our method with the standard SUPG virtual element discretization ^[20] in terms of relative energy errors. In the third test we consider a classic problem taken from ^[31] involving approximation of internal and boundary layers. In all the numerical tests our aim is to assess the robustness of the approach in case of advection dominated regime. Therefore, all the benchmark problems are characterized by large mesh Péclet numbers. Moreover, for each type of polygon in each mesh considered in the tests, we have done a preliminary assessment of the minimum $\ell_E$ that satisfies Assumption 4.1, as detailed in the following section. We consider SFVEM of different orders from one to four and a unit square domain $\Omega = (0, 1) \times (0, 1)$.

6.1.   Test 1

In this first test, we consider an advection-diffusion problem characterized by a diffusivity parameter $\varepsilon = 10^{-9}$ and a transport velocity field $\beta = (1, 0.545)$. The size of the meshes is chosen to guarantee that for the selected value of $\varepsilon$ the mesh Péclet number is much greater than one for all $k$. The forcing term $f$ and the boundary conditions are such that the exact solution (depicted in Figure 1) is

where $c_1 = \frac{3}{\sqrt{2 \pi}}$, $c_2 = \frac{1}{2}$, $c_3 = 1000$ and $c_4 = \frac{1}{3.3} \cdot 10^3$. In Figure 1, we can see that the solution exhibits a strong boundary layer in a direction approximately perpendicular to the direction of the transport velocity field $\beta$.

We consider three different families of meshes ($\mathcal{T}_{1}$, $\mathcal{T}_2$ and $\mathcal{T}_3$) and four different refinements for each one of them. The first mesh of each sequence is reported in Figure 2.

The mesh family $\mathcal{T}_1$ consists of standard cartesian elements, the mesh family $\mathcal{T}_2$ is composed of both concave and convex polygons and the mesh family $\mathcal{T}_3$ is created by Polymesher ^[35]. The first two groups of meshes are refined splitting the existing elements in half, while the mesh family $\mathcal{T}_3$ is refined by Polymesher. Consequently, the tessellations of this family include polygons having different numbers of edges. In Table 1 we report the mean mesh Péclet number for the first mesh and the last mesh of each mesh family.

Regarding the choice of $\ell_E$, we follow the approach in ^[27,28]. On each polygon $E$ of each mesh we perform an eigenvalue analysis on the local matrix

and select the smallest value of $\ell_E$ such that the numerical rank of the matrix is $N_{\mathrm{dof}, k}^E-1$, $N_{\mathrm{dof}, k}^E$ being the number of local degrees of freedom. In particular, we select $\ell_E$ such that the matrix $A^E$ has $N_{\mathrm{dof}, k}^E-1$ eigenvalues that are greater than $1\mathrm{e}-8$. The resulting values are shown in Table 2. Notice that $\ell_E$ does not depend only on the number of vertices of the polygon, but also on its geometry. Indeed, if we consider the line corresponding to $k = 2$ in Table 2, we can see that we require $\ell_E = 2$ for $\mathcal{T}_1$, where quadrilaterals are all squares, and $\ell_E = 1$ for $\mathcal{T}_3$, that features generally shaped quadrilaterals.

In Figures 3–5 we show the convergence curves in log-log scale for orders from $1$ to $4$. We report the relative errors computed in the energy norm plotted against the maximum diameter of the discretization for both the SFVEM and the classical VEM ^[20]. The computed relative error is based on the difference between the exact solution and the projection $\Pi_k^{\nabla, E}$ of the discrete solution $u_{h}$ and it is given by the following expression

The plots in Figures 3–5 display two y-axis: the one on the left is related to the relative energy error, whereas the one on the right is related to the ratio between the VEM error and the SFVEM error. In the legends, alpha denotes the numerical rates of convergence, computed using the last two meshes of each refinement. The numerical rates of convergence for both the two methods are in agreement with the theoretical findings for the energy norm of the problem (Theorem 5.2).

Figure 3a shows that for order $k = 1$ the two methods provide very close results on the mesh family $\mathcal{T}_1$, whereas Figures 4a and 5a show that the SFVEM performs better than the classical VEM on the mesh families $\mathcal{T}_2$ and $\mathcal{T}_3$. Thus, the results suggest that the SFVEM is able to decrease the magnitude of the error with respect to the classical VEM when dealing with solutions characterized by strong anisotropies. An analogous trend is observed for $k = 2$. Moreover, in Figure 4a, 4b, we notice that the error difference between the SFVEM and the VEM is stronger than the one observed in Figures 3a, 3b and 5a, 5b.

Figures 3c, 3d, 4c, 4d, and 5c, 5d show that for $k = 3$ and $k = 4$ the VEM and the SFVEM exhibit an almost equivalent behaviour for all the considered meshes. A similar trend was also observed for anisotropic elliptic problems in ^[24]. Indeed, for higher orders we expect that the polynomial part of the standard VEM bilinear form is more relevant than the stabilizing part, reducing the error gap between the two methods.

6.2.   Test 2

In this second test, we consider an advection-diffusion problem characterized by a diffusivity parameter $\varepsilon = 10^{-9}$ and a transport velocity field $\beta = (1, 0)$. The forcing term $f$ and the boundary conditions are such that the exact solution is $u(x, y) = \sin(2\pi x)\sin(80\pi y)$. We compare the relative energy errors of VEM and SFVEM only on $\mathcal{T}_3$. Indeed, the results of the two methods coincide on $\mathcal{T}_1$ and $\mathcal{T}_2$ since they are by construction aligned with the principal directions of the error (see for details ^[36]). Concerning the choice of $\ell_E$ for the application of SFVEM, we follow Table 2.

In Figure 6, we report the convergence plots for both methods. The results in the lowest order case show a larger difference between the two methods with respect to the previous test. As observed in ^[24], this is due to the highly oscillating behaviour of the solution.

6.3.   Test 3

For the third test, we consider a classic benchmark problem in the SUPG literature characterized by the presence of different layers. This problem was originally proposed in the context of Finite Element Methods in ^[31]. The computational domain $\Omega$ as well as the boundary conditions are depicted in Figure 7. Notice the discontinuous Dirichlet boundary condition on the left side of the domain. We set $\varepsilon = 10^{-6}$ and $\beta = (\cos{\theta}, \sin{\theta})$, where $\theta = \frac{\pi}{4}$, and the forcing term $f$ is zero. The solution features an internal boundary layer due to the discontinuity of the Dirichlet boundary conditions and the highly convective regime, whereas at the outflow boundary another steep layer is produced due to the homogeneous boundary conditions.

We solve the problem using tessellation $\mathcal{T}_2$, depicted in Figure 2b. The resulting Péclet number is very large (around $10^6$). In Figure 8 we display the discrete solutions obtained by the SFVEM scheme and the standard VEM scheme. As it is typical for this problem we notice the presence of undershoots and overshoots near the internal boundary layer.

7.   Conclusions

We presented a numerical investigation of the performances of SUPG-stabilized stabilization-free VEM, assessing the possibility of using higher-order polynomial projectors in the definition of the discrete bilinear form and avoiding the use of a non-polynomial stabilizing bilinear form that does not preserve the structure of the problem operator. We also provided an interpolation estimate for the new scheme, analogous to the one obtained for standard VEM. Numerical results show that the possibility of avoiding a non-polynomial stabilizing bilinear form can enhance the performances of low order classical VEM methods in the case of advection-dominated problems. Indeed, we observed a reduced magnitude of the error especially comparing the lowest order methods, while similar behaviours are observed when choosing higher orders.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors are members of the INdAM-GNCS. A. Borio and F. Marcon kindly acknowledge financial support provided by the European Union through PNRR M4C2 project of CN00000013 National Centre for HPC, Big Data and Quantum Computing (HPC) CUP:E13C22000990001. A. Borio kindly acknowledges partial financial support provided by INdAM-GNCS Projects 2022 and 2023 and by the PRIN 2020 project (No. 20204LN5N5_003). A. Borio also acknowledges financial support by the European Union through project Next Generation EU, M4C2, PRIN 2022 PNRR project CUP:E53D23017950001.

Conflict of interest

The authors declare no conflicts of interest.