The mixed virtual element discretization for highly-anisotropic problems: the role of the boundary degrees of freedom

1.   Introduction

The Virtual Element Method (in short VEM) ^[6,7] is a generalization of the Finite Element Method (FEM in short) that can easily handle general polytopal meshes and high-order methods. The major difference with the FEM is that the VEM space contains suitable non-polynomial functions. For this reason, the standard VEM discrete bilinear form is the sum of a consistency part ensuring accuracy and a stabilization term enforcing the coercivity. In particular, the choice of the stabilization term remains a critical part of the VEM construction ^[12,26] and it is usually problem-driven. Furthermore, the stabilization term may have possible negative effects on the conditioning of the system ^[16,24] and may become an issue in highly anisotropic diffusion problems due to its isotropic nature. Moreover, we recall that it is crucial to employ a well-conditioned polynomial basis in the definition of the internal degrees of freedom to obtain reliable solutions when building high-order methods. Indeed, the advantages of using $L^2$-orthogonal polynomial bases against the standard monomial one have largely been proved both for the primal version of the method ^{[5,11,13,16,24]} and for its mixed formulation ^[14].

In this paper, we want to test the accuracy and the robustness of the mixed virtual element Method when dealing with highly anisotropic diffusion tensors. For this purpose, we propose different kinds of degrees of freedom and test them against different choices of the stabilization term for a set of benchmark anisotropic diffusion problems. In particular, we introduce a new set of boundary degrees of freedom which are defined as moments up to degree $k \geq 0$ against an $L^2([0, 1])$-orthonormal polynomial basis in order to analyze the role of the boundary degrees of freedom in the conditioning and in the accuracy of the methods. Numerical experiments show that this choice of boundary degrees of freedom generally leads to a downward shift of the error curves. However, this approach does not ensure an improvement of the condition number of the system matrix in all the tested cases.

The outline of the paper is as follows. In Section 2 we present the model problem. In Section 3, after introducing the local mixed virtual element spaces and different sets of the local degrees of freedom, we define the mixed VE formulation of the problem. In Section 4, we describe the main properties and discuss possible choices for the stabilization term. Finally, in Section 5 we test all the proposed approaches through different benchmark problems which are characterized by highly anisotropic diffusion tensors, with both constant and variable coefficients.

2.   The model problem

Let $\Omega \subset \mathbb{R}^2$ be a bounded convex polytopal domain with boundary $\Gamma$ and let $\mathit{\boldsymbol{n}}_{\Gamma}$ be the outward unit normal vector to the boundary. Let us consider a tensor $\mathit{\boldsymbol{D}}(\mathit{\boldsymbol{x}}) \in \mathbb{R}^{2\times 2}$ which is bounded, measurable, symmetric and strongly elliptic on $\Omega$, i.e., there exist $\mathit{\boldsymbol{D}}_{min}, \ \mathit{\boldsymbol{D}}_{max}$, independent on $\mathit{\boldsymbol{v}}$ and $\mathit{\boldsymbol{x}}$, such that

holds for every $\mathit{\boldsymbol{v}} \in H_{0, \Gamma_N}(\mathrm{div}; \Omega) = \{v \in H(\mathrm{div}; \Omega): \mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{n}}_{\Gamma_N} = 0\}$ and for almost every $\mathit{\boldsymbol{x}} \in \Omega$, where ${\lVert \cdot \rVert}$ denotes the euclidean norm. Given $f\in L^2(\Omega)$, $g_D \in H^{1/2}(\Gamma_D)$ and $g_N \in L^2(\Gamma_N)$, we consider the following diffusion problem

where $\Gamma_D$ and $\Gamma_N$ such that $\Gamma_D \cup \Gamma_N = \Gamma$, $\vert \Gamma_D \vert \neq 0$ and $\vert \Gamma_D \cap \Gamma_N \vert = 0$ denote the Dirichlet and the Neumann boundary, respectively. In particular, in the following, we focus on diffusion problems with a diffusion tensor of the form

which is characterized by a high anisotropic ratio, i.e., the ratio between the smallest and largest eigenvalues of the diffusion tensor.

Introducing the velocity space $\mathit{\boldsymbol{V}} = H_{0, \Gamma_N}(\mathrm{div}; \Omega)$ and the pressure space $Q = L^2(\Omega)$, the mixed variational formulation of (2.1) reads:

where $\mathit{\boldsymbol{u}}_N \in H(\mathrm{div}; \Omega)$ is a chosen function that satisfies $\mathit{\boldsymbol{u}}_N \cdot \mathit{\boldsymbol{n}}_{\Gamma_N} = g_N$ and $\langle \cdot, \cdot \rangle_{\pm\frac{1}{2}, \Gamma_D}$ denotes the duality paring between $H^{-1/2}(\Gamma_D)$ and $H^{1/2}(\Gamma_D)$.

3.   The mixed virtual element space

Now, let us consider a decomposition $\mathcal{T}_{h}$ of $\Omega$ in star-shaped polygons $E$, where $h$, as usual, is set to be the maximum diameter of elements $E \in \mathcal{T}_{h}$. We further denote by $\mathcal{E}_{h, E}$ the set of edges of an element $E \in \mathcal{T}_{h}$.

For any integer $k \geq 0$, we define the local virtual element space related to the velocity variable $\mathit{\boldsymbol{u}}$ as

and the local virtual element space related to the pressure variable $p$ as $\mathbb{P}^{}_{{k}}\!\left({E}\right)$.

The choice of the degrees of freedom in the local pressure space $\mathbb{P}^{}_{{k}}\!\left({E}\right)$ is trivial: The degrees of freedom of a function $p \in \mathbb{P}^{}_{{k}}\!\left({E}\right)$ are its coefficients with respect to the polynomial basis chosen as the basis for $\mathbb{P}^{}_{{k}}\!\left({E}\right)$. The standard polynomial basis for $\mathbb{P}^{}_{{k}}\!\left({E}\right)$ used in the VEM construction ^[9] is given by the set of the $n_{k} = \dim \mathbb{P}^{}_{{k}}\!\left({E}\right) = \frac{(k+1)(k+2)}{2}$ bi-dimensional scaled monomials, i.e.,

where $\mathit{\boldsymbol{x}}_E$ and $h_E$ are the centroid and the diameter of the polygon $E$, respectively, and $\ell: \mathbb{N}^2 \to \mathbb{N}$ is the function that maps

A more robust choice is represented by the set of the $L^2(E)$-orthonormal polynomials $\mathcal{Q}_{{k}}\!\left({E}\right) = \{q_{\alpha}\}_{\alpha = 1}^{n_k}$ introduced in ^[5,11,24] for the primal version of the method and then tested in the mixed case in ^[14]. This orthonormal polynomial basis is defined as

where ${\bf L}^k \in \mathbb{R}^{n_k \times n_k}$ is built by applying twice the modified Gram-Schmidt algorithm to the monomial Vandermonde matrix related to a proper quadrature formula on $E$.

3.1.   The degrees of freedom for the velocity variable

In order to define the local degrees of freedom for the local velocity space $\mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}$, we need to introduce the following polynomial spaces. We introduce the (vector) polynomial space

and the set $\mathcal{G}_{{k}}^{{\perp, m}}\!\left({E}\right) = \Big\{ \mathit{\boldsymbol{g}}^{\perp, m}_{\alpha}\Big\}_{\alpha = 1}^{n^{\perp}_k}$ which is defined in such a way

with $n^{\nabla}_k = n_k + (k+1)$ and $n^{\perp}_k = n_{k}- (k+1)$. The set

represents a (vector) polynomial basis for $\left[\mathbb{P}^{}_{{k}}\!\left({E}\right)\right]^2$ which allows to easily define the set of local degrees of freedom in the mixed VEM framework ^[9]. Let us denote by

the set of (vector) polynomials

such that

which is obtained by orthonormalizing the gradients of polynomials belonging to $\mathcal{Q}_{{k+1}}\!\left({E}\right)$ through the modified Gram-Schmidt algorithm. Now, we define $\mathcal{G}_{{k}}^{{\perp, \overline{q}}}\!\left({E}\right) = \{ \mathit{\boldsymbol{g}}^{\perp, \overline{q}}_{\alpha}\}_{\alpha = 1}^{n^{\perp}_{k}}$ as the $L^2(E)$-orthogonal complement of $\mathcal{G}_{{k}}^{{\nabla, \overline{q}}}\!\left({E}\right)$ in $\left[\mathbb{P}^{}_{{k}}\!\left({E}\right)\right]^2$, which is chosen such that

Further details about the construction of this basis can be found in ^[14]. Here, it was shown that it is advisable to choose the set

as the (vector) polynomial basis for $\left[\mathbb{P}^{}_{{k}}\!\left({E}\right)\right]^2$ in order to reduce the ill-conditioning of the system matrix and to obtain more accurate and reliable solutions for high values of the local polynomial degree and in the presence of badly-shaped polygons.

Now, let us introduce a quadrature formula $\mathbb{S}^Q = \{(s_j^Q, w_j^Q)\}_{j = 1}^{N^Q}$ of order $2(k+1)$ with $N^Q \geq k+2$ nodes on the interval $[0, 1]$. We define the one-dimensional $L^2([0, 1])$-orthonormal polynomial basis $\mathcal{Q}_{{k+1}}\!\left({[0, 1]}\right) = \{t_1, \dots, t_{k+1}, t_{k+2}\}$ for $\mathbb{P}^{}_{{k+1}}\!\left({[0, 1]}\right)$ by applying the modified Gram-Schmidt algorithm with reorthogonalization to the Vandermonde matrix ${\bf V}^{\mathbb{S}^Q} \in \mathbb{R}^{N^Q \times (k+2)}$ related to the one-dimensional monomial basis $\{1, s, \dots, s^k, s^{k+1}\}$ and the quadrature formula $\mathbb{S}^Q$. More precisely, we perform sequentially

where ${\bf W}^{\mathbb{S}^Q} \in \mathbb{R}^{N^Q \times N^Q}$ is the diagonal matrix of quadrature weights, and then we define

where ${\bf L}^{\mathbb{S}^Q} = ({\bf R}^{\mathbb{S}^Q}_2 {\bf R}^{\mathbb{S}^Q}_1)^{-T}$.

We remark that each polynomial in $\mathbb{P}^{}_{{k+1}}\!\left({e}\right)$, $e\in\mathcal{E}_{h, E}$, can be written in terms of polynomials in $\mathcal{Q}_{{k+1}}\!\left({[0, 1]}\right)$ through an affine mapping $F_e:[0, 1] \to e$. Furthermore, we recall that the modified Gram-Schmidt algorithm is a hierarchical procedure, which means, for example, that

is a basis for $\mathbb{P}^{}_{{k}}\!\left({[0, 1]}\right)$.

In $\mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}$, we define the set of local Degrees of Freedom (DOFs in short) as the union of

(1) the set of the boundary degrees of freedom which can be chosen as

(1)a the values of $\mathit{\boldsymbol{v}}_h \cdot \mathit{\boldsymbol{n}}_e$ in the $k+1$ Gauss quadrature points $\mathit{\boldsymbol{x}}_i^{e, Q}$ internal on each edge $e \in \mathcal{E}_{h{, E}}$,

or

(1)b the $k+1$ moments on each edge $e \in \mathcal{E}_{h{, E}}$:

where $\vert e \vert$ represents the length of the edge $e$;

(2) the set of the internal degrees of freedom which can be chosen as the internal moments computed against

(2)a the sets of functions $\mathcal{G}_{{k-1}}^{{\nabla, m}}\!\left({E}\right)$ and $\mathcal{G}_{{k}}^{{\perp, m}}\!\left({E}\right)$:

or

(2)b the sets of functions $\mathcal{G}_{{k-1}}^{{\nabla, \overline{q}}}\!\left({E}\right)$ and $\mathcal{G}_{{k}}^{{\perp, \overline{q}}}\!\left({E}\right)$:

where $\vert E \vert$ is the area of the polygon $E$.

Let us denote by $N^{\operatorname{dof}}_E = \dim \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)} = \#\mathcal{E}_{h{, E}}(k+1)+ n^{\nabla}_{k-1} + n^{\perp}_k$ and let us introduce the local Lagrangian VE basis $\{\mathit{\boldsymbol{\varphi}}_i\}_{i = 1}^{N^{\operatorname{dof}}_E}$ related to the local degrees of freedom, where the DOF numbering first counts the boundary DOFs and then the internal DOFs. Furthermore, for each element $E \in \mathcal{T}_{h}$, we define the operators $\operatorname{dof}_i : \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}\to \mathbb{R}$ which associate each function $\mathit{\boldsymbol{v}} \in \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}$ to its $i$-th degree of freedom.

Now, let us introduce the $L^2(E)$-projector $\mathit{\boldsymbol{\Pi}}^{{0, E}}_{{k}}: \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}\to \left[\mathbb{P}^{}_{{k}}\!\left({E}\right)\right]^2$, which is defined by the orthogonality condition

We note that each combination of the aforementioned degrees of freedom makes the projection $\mathit{\boldsymbol{\Pi}}^{{0, E}}_{{k}} \mathit{\boldsymbol{v}}_h$ of a function $\mathit{\boldsymbol{v}}_h \in \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}$ computable. In particular, the computation of $\mathit{\boldsymbol{\Pi}}^{{0, E}}_{{k}} \mathit{\boldsymbol{v}}_h$ with the pairs consisting of DOFs (1)a and (2)a, and DOFs (1)a and (2)b has been largely discussed in ^[9,14]. Concerning the choice of DOFs (1)b and (2)a, we first note that, given $\mathit{\boldsymbol{v}} \in \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}$, the orthogonality condition (3.12) yields

and

The right-hand side of Eq (3.14) can be computed through the internal degrees of freedom (3.9). Now, we recall that $\mathrm{div} \mathit{\boldsymbol{v}}$ is a polynomial $\sum_{\alpha = 1}^{n_k} c_{\alpha} m_{\alpha}^k \in \mathbb{P}^{}_{{k}}\!\left({E}\right)$ whose coefficients $\{c_{\alpha}\}_{\alpha = 1}^{n_k}$ can be determined by imposing

The first term of the right-hand side of (3.15) can be computed through the internal degrees of freedom (3.8). Furthermore, on each edge $e \in \mathcal{E}_{h, E}$, we can write the trace of monomials as

Thus, the second term of the right-hand side of (3.15) becomes

which can be computed by resorting to the boundary degrees of freedom (1)b.

In order to compute the second term of the right-hand side of Eq (3.13), we should determine the polynomial $\mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{n}}_e$ on each edge $e \in \mathcal{E}_{h, E}$. Nonetheless, if $\{\mathit{\boldsymbol{\varphi}}_i^e\}_{i = 1}^{k+1}$ is the local Lagrangian mixed VE basis related to the boundary degrees of freedom defined on the edge $e\in \mathcal{E}_{h, E}$, we observe that

while $\mathit{\boldsymbol{\varphi}}\cdot \mathit{\boldsymbol{n}}_e$ is the zero-polynomial on $e$ if $\mathit{\boldsymbol{\varphi}}$ is related to an internal degree of freedom or to a different edge of $E$. Finally, since $\mathcal{Q}_{{k+1}}\!\left({[0, 1]}\right)$ is an $L^2([0, 1])$-orthonormal basis for $\mathbb{P}^{}_{{k+1}}\!\left({[0, 1]}\right)$, we simply have $\forall \alpha = 1, \dots, n_k^{\nabla}$, $i = 1, \dots, k+1$ and $\forall e \in \mathcal{E}_{h, E}$

The construction of the method with the choice of DOFs (1)b and (2)b is analogous to the one which exploits the degrees of freedom (1)b and (2)a. Indeed, we recall that we are able to write

where ${\bf L}^{\nabla, k}$ and ${\bf L}^{k+1}$ are defined in (3.4) and (3.2), respectively.

Remark 3.1. Note that, since we define the one-dimensional polynomial basis $\mathcal{Q}_{{k}}\!\left({[0, 1]}\right)$ on the interval $[0, 1]$, we must perform the orthogonalization process just once. Thus, the additional cost in taking an $L^2([0, 1])$-orthonormal basis instead of the one-dimensional monomial basis is negligible and independent of the number of edges of the tessellation $\mathcal{T}_{h}$.

3.2.   The mixed virtual element formulation of the model problem

On each element $E \in \mathcal{T}_{h}$, let us define the continuous local bilinear form

and its discrete counterpart

which is the sum of the consistency term

and of the stability term $S^{{{E}}} {\left(\cdot, \cdot\right)}$, which is any symmetric positive-definite bilinear form that satisfies

for some positive constants $\alpha_{\ast}, \ \alpha^{\ast}$ depending on $\mathit{\boldsymbol{D}}^{-1}$ but independent on $h$ ^[8,15].

Now, let us introduce the global mixed virtual element spaces

for the velocity and the pressure variables, respectively. In particular, as global degrees of freedom for each $\mathit{\boldsymbol{v}}_h \in \mathit{\boldsymbol{V}}_{h, {k}}$, we consider

● the boundary degrees of freedom of $\mathit{\boldsymbol{v}}_h$ defined on each internal edge of the decomposition and at edge boundary with Dirichlet boundary conditions;

● the internal degrees of freedom in each element $E\in \mathcal{T}_{h}$.

Furthermore, the value of the boundary DOFs at the Neumann edge is fixed in accordance with the value of the Neumann boundary conditions.

Finally, the virtual element discretization of the problem (2.3) reads

where $\mathit{\boldsymbol{u}}_{N, h} \in \Big\{ \mathit{\boldsymbol{v}} \in H(\mathrm{div}; \Omega): \mathit{\boldsymbol{v}} \in \mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}\forall E \in \mathcal{T}_{h}\Big\}$ is such that

for each boundary degree of freedom $i$ related to a Neumann edge.

4.   The stabilization term

Let us introduce the elemental matrix ${\bf A}^{E} \in \mathbb{R}^{N^{\operatorname{dof}}_E\times N^{\operatorname{dof}}_E}$, whose entries are defined as the application of the local discrete bilinear form $a_h^{{{E}}} {\left(\cdot, \cdot\right)}$ to the Lagrangian basis functions of $\mathit{\boldsymbol{V}}_{h, {k}} {\left(E\right)}$, i.e., $\forall i, j = 1, \dots, N^{\operatorname{dof}}_E$

where ${\bf A}^{E}_C$ and ${\bf A}^{E}_S$ represent the elemental matrices related to the consistency and the stability term, respectively. The complete elemental matrix related to the mixed discretization of the problem (3.22) reads

where the entries of the divergence matrix ${\bf W}^E \in \mathbb{R}^{n_k \times N^{\operatorname{dof}}_E}$ are defined as

Since the degrees of freedom of the velocity space are chosen in such a way the related Lagrangian VE basis functions scale uniformly with respect to the mesh size $h$, the most natural mixed VEM stabilization $S^{{{E}}} {\left(\cdot, \cdot\right)}$ which satisfies (3.21) is the so-called dofi-dofi stabilization ^[7,8]:

where $C_{ \mathit{\boldsymbol{D}}^{-1}}$ is a constant depending on $\mathit{\boldsymbol{D}}^{-1}$. Moreover, $\forall \mathit{\boldsymbol{u}} \in \mathit{\boldsymbol{V}}_{h, {k}}$, from the definition of the internal degrees of freedom and the definition (3.12) of the projector $\mathit{\boldsymbol{\Pi}}^{{0, E}}_{{k}}$, it follows

for each internal degree of freedom $i$. Thus, in the mixed VEM construction, it is not necessary to include the internal degrees of freedom in the stabilization procedure ^[8].

Furthermore, as highlighted in ^[10], in order to avoid to level off the stabilization term with respect to the consistency term for the higher polynomial degrees, which would lead to a loss of accuracy, we can choose the so-called D-recipe stabilization, defined as follows

where $S_{ii} = C_{ \mathit{\boldsymbol{D}}^{-1}} \vert E \vert \max(1, ({\bf A}^E_C)_{ii})$ if $i$ is related to a boundary degree of freedom and $S_{ii} = 0$ otherwise, since we do not need to include the internal degrees of freedom (see Eq (4.2)).

Usually, the constant $C_{ \mathit{\boldsymbol{D}}^{-1}}$ is taken equal to the spectral norm $\| \mathit{\boldsymbol{D}}^{-1}\| = 1/ \mathit{\boldsymbol{D}}_{min}$, since $\mathit{\boldsymbol{D}}$ is assumed to be symmetric and strongly elliptic.

Finally, the choice of the stabilization term and, in particular, of the constant $C_{ \mathit{\boldsymbol{D}}^{-1}}$ should be dependent on the problem features and on the definition of the local degrees of freedom ^[7,26].

5.   Numerical experiments

In this section, we perform some numerical experiments that allow us to show the role of the boundary degrees of freedom and of the stabilization term in preventing the ill-conditioning of the system matrix. To this end, we analyze the behaviour of the condition number of the global system matrix ${\bf K}$ and of the errors

at varying of the polynomial degree $k$ or of the mesh size $h$, for different families of meshes. Given $k \geq 0$ and the mesh size $h$, we recall that if the solution is sufficiently smooth, the expected convergence rates of errors (5.1) and (5.2) are $O(h^{k+1})$.

In the following, we use the notation:

● Mon (a) to denote the approach which exploits the pair of DOFs (1)a and (2)a;

● Mon (b) to denote the approach which exploits the pair of DOFs (1)b and (2)a;

● Ortho (a) to denote the approach which exploits the pair of DOFs (1)a and (2)b;

● Ortho (b) to denote the approach which exploits the pair of DOFs (1)b and (2)b.

We note that we use the scaled monomial basis as the basis for the local pressure space in the monomial approaches (Mon), while we use the $\mathcal{Q}_{{k}}\!\left({E}\right)$ basis in the orthonormal approaches (Ortho).

5.1.   Test 1: boundary degrees of freedom

In this first test, we analyze the behaviour of the four aforementioned approaches by solving a Poisson problem with homogeneous Dirichlet boundary conditions.

More precisely, let us set $\Omega = (0, 2)^2$, we define the forcing term $f$ in such a way the exact pressure is

In this test, we employ the dofi-dofi stabilization term with $C_{ \mathit{\boldsymbol{D}}^{-1}} = 1$ and we evaluate the performance of our approaches using a set of three concave meshes $\{\mathcal{T}_{h{_{i}}}^C\}_{i = 1}^3$. These meshes are generated through an agglomeration process that begins with basic convex meshes at different refinement levels. The second refinement $\mathcal{T}_{h{_{2}}}^C$ is depicted in Figure 1.

In Figure 2, we show the behaviour of the condition number of the global system matrix ${\bf K}$ at varying of the polynomial degree $k$, for each concave mesh $\mathcal{T}_{h{_{i}}}^C$, $i = 1, 2, 3$, in semilog plots. From these graphs, we note that changing the boundary degrees of freedom from (1)a to (1)b generally does not ensure an improvement in the condition number of the global system matrix for fixed internal degrees of freedom. Furthermore, we observe that, in order to cure the ill-conditioning of the global system matrix, the use of an $L^2(E)$-orthonormal (vector) polynomial basis for $\left[\mathbb{P}^{}_{{k}}\!\left({E}\right)\right]^2$ is strongly recommended, as already highlighted in ^[14].

Figures 3 and 4 show the behaviour of errors (5.1) and (5.2) at varying of the polynomial degree $k$ for each $\mathcal{T}_{h{_{i}}}^C$, with $i = 1, 2, 3$, in semilog plots. Furthermore, Figures 5 and 6 show the behaviour of such errors for decreasing values of the mesh size $h_i$, $i = 1, 2, 3$, for $k = 1, 3, 5$, with a loglog scale. From these figures, we can note that changing the internal degrees of freedom from (2)a to (2)b does not significantly modify the behaviour of errors (5.1) and (5.2), at varying of the mesh size $h$, for the lower values of the polynomial degree $k$. In general, this is not true for the boundary degrees of freedom. Indeed, from Figures 3 and 5, we can note that the error (5.1) is sensitive to a variation from (1)a to (1)b of the boundary degrees of freedom, especially on the coarser meshes. As the mesh is refined, such difference becomes smaller and smaller and the orthonormal approaches tend to behave in the same way regardless of the type of boundary DOFs used.

Finally, for the higher values of $k$, the errors start to raise due to the ill-conditioning of the matrix ${\bf K}$ in the Mon approaches, while the Ortho approaches are robust also for the higher polynomial degrees.

5.2.   Test 2: anisotropic diffusion problems

In this experiment, we want to analyze the sensitivity of the presented approaches to the choice of stabilization in a context where such sensitivity becomes the main issue to overcome, namely the diffusion problems with high anisotropic coefficients.

Equations characterized by anisotropic diffusion coefficients arise in many practical contexts, such as the heat equation, groundwater flow, transport problems and so on. Generally, these types of problems are expressed as parametric problems and they are numerically treated by means of ad hoc methods, needed to avoid the so-called locking phenomenon ^[4]. This phenomenon occurs experimentally when the discretization error does not decrease at the expected rate when the parameter tends to limiting values and, in general, it is typical of the lower order schemes. These ad hoc methods include variational crimes, i.e., modification of the bilinear form ^[19], and flow-aligned grid methods ^[21]. In particular, in the Virtual Element context, the isotropic nature of the standard stabilization term can become an issue in these kinds of problems and different approaches have been studied to handle the anisotropic nature of the diffusion tensors ^[12,25] mainly for the primal formulation of the method.

Thus, we consider the test problem proposed in ^[23], which is a dimensionless parametric version of problem (2.1) with a constant diffusion tensor, defined on $\Omega = (0, 1)^2$. In particular, the diffusion tensor $\mathit{\boldsymbol{D}} = \begin{bmatrix} 1 & 0 \\ 0 & \epsilon \end{bmatrix}$ depends on the diffusion parameter $\epsilon \in [10^{-6}, 1]$, which, in this case, represents also the anisotropic ratio. In our notation, $\mathit{\boldsymbol{D}}_{min} = \epsilon$ (or $\mathit{\boldsymbol{D}}^{-1}_{max} = \frac{1}{\epsilon}$) and $\mathit{\boldsymbol{D}}_{max} = 1$.

The performance of the four approaches is evaluated on two different kinds of families of meshes: a cartesian $\mathcal{T}_{h}^Q$ family and a family $\mathcal{T}_{h}^{DQ}$ of distorted quadrilateral meshes generated from the cartesian ones through a sine distortion. For each family of meshes, we consider four refinements $\{\mathcal{T}_{h{_{i}}}^Q\}_{i = 1}^{4}$ and $\{\mathcal{T}_{h{_{i}}}^{DQ}\}_{i = 1}^{4}$: the first and the last refinement of each family are shown in Figure 7.

In order to compute errors (5.1) and (5.2), we choose the parametric exact solution

The presence of $\epsilon$ at the exponent of (5.3) makes the low conductivity direction dominant when $\epsilon$ tends to zero and the nearly pure Neumann boundary conditions are set, by leading, in general, to very poor results when employing standard methods ^[19]. Thus, we test three different kinds of boundary conditions (BCs in short):

● pure Dirichlet boundary conditions, i.e., $\Gamma_D = \Gamma$;

● mixed Dirichlet-Neumann boundary conditions with $\Gamma_D = \{(x, y): x = 0 \text{ or } y = 0\};$

● nearly pure Neumann conditions, that is we set $\Gamma_D = \{(x, y): (x = 1 \text{ and } 1-\delta \leq y \leq 1) \text{ or } (y = 1 \text{ and } 1-\delta \leq x \leq 1)\},$ where $\delta$ decreases with the mesh size as $\frac{1}{5 \cdot 2^{i-1}}$ $i = 1, \dots, 4$.

In the first two cases, generally, no locking phenomenon occurs.

Furthermore, we test three possible choices for the stabilization term, namely

● S$_1$: the standard dofi-dofi stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = \| \mathit{\boldsymbol{D}}^{-1} \| = \frac{1}{\epsilon}$;

● S$_2$: the standard dofi-dofi stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = 1$;

● S$_3$: the D-recipe stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = 1$.

We observe that when $\epsilon$ becomes very small, the constant $C_{ \mathit{\boldsymbol{D}}^{-1}}$ related to the choice S$_1$ becomes very big.

5.2.1.   Effect of the anisotropy on the condition number of the global system matrix

In Figures 8 and 9 we report the behaviour of the condition number of the global system matrix at varying of $k$ in semilog plots, when the Dirichlet and the nearly pure Neumann boundary conditions are set, respectively. The results are related to $\epsilon \in \{1, 10^{-6}\}$ and to the $\mathcal{T}_{h{_{1}}}^{Q}$ and the $\mathcal{T}_{h{_{1}}}^{DQ}$ meshes.

Accordingly to results presented in ^[14], we observe an exponential growth in the condition number of the matrix ${\bf K}$ when the internal DOFs (2)a are employed. A linear growth is observed instead when resorting to the choice (2)b. Furthermore, as already pointed out in the previous test, changing the boundary DOFs from (1)a to (1)b does not generally lead to an improvement of the behaviour of the condition number of ${\bf K}$.

Furthermore, we note that a sine distortion of elements causes a faster increase in the condition number of ${\bf K}$ when the internal DOFs (2)a are used, while this growth is not so evident in the case of the internal DOFs (2)b.

We further note that having nearly pure Neumann boundary conditions has just a small effect on the condition number of ${\bf K}$ for the lower values of $k$ and that the condition number of ${\bf K}$ seems to be mainly controlled by the anisotropic effect, as observed in ^[23].

Finally, we observe that the pair of DOFs (1)b and (2)b reveals to be the most robust approach with respect to the choice of the stabilization term, whereas stabilization choice S$_1$ seems to be the worst choice in terms of the condition number of ${\bf K}$, if a combination of DOFs different from (1)b and (2)b is used.

5.2.2.   The mesh alignment and the locking phenomenon

Figures 10–12 show the behaviour of the pressure error (5.1) at varying of the polynomial degree $k$ related to the Dirichlet, mixed and nearly pure Neumann boundary conditions, respectively. These results are obtained on the cartesian mesh $\mathcal{T}_{h{_{1}}}^{Q}$.

From these figures we observe that, after an initial decrease, the error starts to raise due to ill-conditioning, but only when the internal DOFs (2)a are employed. Choosing internal DOFs (2)b leads to the best performance in each tested case for the higher values of the polynomial degree $k$.

The error curves related to the different analyzed approaches are very similar for the lower values of $k$ when a cartesian mesh is used. The only exception is represented by the choice of the boundary DOFs (1)a and stabilization term S$_1$. In this case, error curves are slightly upward shifted for the smaller values of $\epsilon$ when nearly pure Neumann boundary conditions are set.

In Figures 13–15 we report the behaviour of the pressure error (5.1) at varying of $k$ related to the Dirichlet, mixed and nearly pure Neumann boundary conditions, in the case of the distorted cartesian mesh $\mathcal{T}_{h{_{1}}}^{DQ}$. By comparing these results with those obtained in the case of cartesian mesh, we can observe that, in the case of distorted meshes, the considered approaches show very different behaviours in terms of error (5.1) when $\epsilon$ is very small. Indeed, we highlight that the Cartesian mesh is aligned with the directions of the anisotropy, by limiting the effect of anisotropy. The main variations are observed for the approaches that exploit the stabilization term S$_1$ also in the case of distorted meshes. Furthermore, we must observe an initial upward shift of the error curves related to the D-recipe S$_3$ for the lower values of the polynomial degree $k$ with respect to the approaches that use the stabilization term S$_2$. Nevertheless, for the higher values of $k$, the stabilization terms S$_2$ and S$_3$ yield again similar results and very good performance is achieved when they are combined with the internal DOFs (2)b.

In order to analyze better such differences, in Figures 16 and 17 we report the behaviour of the errors (5.1) and (5.2) at decreasing values of the mesh size $h$ for the lowest polynomial degree $k = 0$ and in the case of pure nearly Neumann conditions for the Cartesian and the distorted quadrilateral families of meshes, respectively. In the lowest-order case, we can observe a locking phenomenon in the pressure error when distorted quadrilateral meshes are employed, as suggested by an upward shift of the error curves when $\epsilon \to 0$ and by a loss in the convergence rates, which can describe a pre-asymptotic regime ^[4,23]. As mentioned before, the locking phenomenon is typical, in general, of the lower order methods. Indeed, looking at Figures 18 and 19, we can note that the approaches which employed orthogonal internal DOFs show the right rates of convergence for the higher values of $k$. The monomial approaches, instead, do not converge due to ill-conditioning when $k$ is high.

5.3.   Test 3: two magnetic islands

In the previous experiment, we considered a constant diffusion tensor with the diffusion directions aligned with the Cartesian axes. Thus, the problem of anisotropy could be easily handled by choosing a Cartesian family of meshes.

Now, we propose the test "Two Magnetic Islands" described in ^[18], where it is almost impossible to generate an aligned mesh to solve the problem. This example models the instability phenomenon which arises in magnetized plasma for fusion applications. More precisely, we consider a diffusion problem in $\Omega = (-1, 1) \times (-0.5, 0.5)$, with a diffusion tensor given by

where the unit vector $\mathit{\boldsymbol{b}} = \begin{bmatrix} b_1 & b_2 \end{bmatrix}^T$ represents the parallel direction to the anisotropy (or to the magnetic field $\mathit{\boldsymbol{B}}$), while $D_{\vert \vert}$ and $D_{\perp}$ represent the parallel and the perpendicular diffusion coefficients, respectively. In this kind of application, we observe that $D_{\vert \vert}$ can be greater than $D_{\perp}$ by a factor of $10^{12}$ ^[18]. Let us now define the equilibrium magnetic field

which is shown in Figure 20(a). By looking at this figure, we note that the magnetic field results to be the zero-vector in the center of the "magnetic islands" (the $\mathcal{O}$-points) and where the field lines cross each other (the so-called $\mathcal{X}$-points). In all the other points, we can define

and compute

We further fix $D_{\perp} = 1$, while $D_{\vert \vert} \in \{1, 10^{4}, 10^{8}\}$.

We evaluate the performance of the aforementioned approaches on a family $\mathcal{T}_{h}^S = \{\mathcal{T}_{h{_{i}}}^S\}_{i = 1}^{4}$ of four squared meshes, which are characterized by an edge length decreasing as $\frac{1}{2^{i+1}}$, with $i = 1, \dots, 4$. We note that both the $\mathcal{O}$-points and $\mathcal{X}$-points represent vertices of the tessellation in each refinement. Furthermore, we define the forcing term and the boundary conditions in such a way the exact solution is

which is shown in Figure 20(b). We test two cases, characterized by different boundary conditions, namely

● pure Dirichlet boundary conditions $\Gamma_D = \Gamma$;

● mixed boundary conditions, with $\Gamma_N = \{(x, y): x = -1 \text{ or } x = 1\}.$

We note that the velocity field does not depend on the parameter $D_{\vert \vert}$.

In this experiment, we test three possible choices for the stabilization term, namely

● S$_1$: the dofi-dofi stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = \frac{1}{D_{\vert \vert}}$.

● S$_2$: the D-recipe stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = 1$.

● S$_3$: a D-recipe stabilization term with

where $\mathit{\boldsymbol{x}}_{e_i}$ and $\mathit{\boldsymbol{n}}_{e_i}$ are the midpoint and the unit outward normal vector to the edge $e_i$ related to the boundary DOF $i$. This stabilization term, inspired by ^[17], aims to take into account the actual strength of the normal contribution of the parallel diffusion on each edge.

Figures 21 and 22 show the behaviour of the errors (5.1) and (5.2) at varying of the polynomial degree $k$ for the second refinement $\mathcal{T}_{h{_{2}}}^S$ when the Dirichlet and mixed boundary conditions are imposed. We decide to report only the behaviour of the Ortho approaches in these figures in order to try to better highlight differences between the employment of boundary DOFs (1)a and (1)b.

We observe that all approaches show the right behaviour in terms of the relative pressure error (5.1) in the case of both Dirichlet and mixed boundary conditions. This appears also evident when observing the behaviour of the pressure error in terms of $h$ in Figures 23–26 for the lowest order $k = 0$ and for the polynomial degree $k = 2$. From these figures we can note that, as usual, the choice of boundary DOFs (1)b is characterized by smaller pressure error constants with respect to the choice (1)a. Furthermore, approaches that employ (1)b seem to be less sensitive to the choice of the stabilization term than approaches which exploit boundary DOFs (1)a.

However, the same conclusions do not hold true when dealing with the relative velocity error (5.2). Indeed, we first can note that switching off the stabilization by choosing the stabilization term S$_1$ when $D_{\vert \vert}$ is big enough generally does not lead to good results in terms of the velocity error. Furthermore, we note that, in order to achieve good results in terms of the velocity error, it is very important to enforce the velocity on the boundary by imposing strong Neumann boundary conditions when high values of $D_{\vert \vert}$ are considered. In this way, it is possible to obtain the right convergence rates in terms of the mesh size for both the pressure and the velocity errors as can be seen in Figures 23 and 24.

Finally, we observe that, in this test case, the Ortho (a) approach seems to perform better than the Ortho (b) approach in terms of the velocity error when highly anisotropic cases are taken into account.

5.4.   Test 4: anisotropic and heterogeneous case over irregular grids

In this last experiment, we propose a benchmark anisotropic diffusion problem ^[20] which is deeply analyzed in ^[3] and inspired by ^[2,22] and which may induce the locking phenomenon on some finite volume schemes in the mixed formulation.

More precisely, let us consider the problem (2.1) on $\Omega = (0, 1)^2$, we define the diffusion tensor as in (5.4), with $D_{\vert \vert} = 1$, $D_{\perp} = 10^{-3}$ and the field

which is represented in Figure 27(a). Thus, in each point different from $(0, 0)$, we have

Let us set the Dirichlet boundary conditions and the forcing term $f$ in such a way the exact solution is $p(x, y) = \sin(\pi x) \sin(\pi y).$

In this test case, we assess the performance of our methods by using six different refinements of an irregular grid, which are available at ^[1], namely, the six refinements of the family of meshes "mesh4_1". The third refinement is shown in Figure 27(b).

In this experiment, we further test four possible choices for the stabilization term, i.e.,

● S$_1$: the dofi-dofi stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = \frac{1}{D_{\perp}}$.

● S$_2$: the dofi-dofi stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = 1$.

● S$_3$: the D-recipe stabilization with $C_{ \mathit{\boldsymbol{D}}^{-1}} = 1$.

● S$_4$: a D-recipe stabilization term with

Figure 28 shows the behaviour of errors (5.1) and (5.2) as the polynomial degree $k$ increases on the coarsest mesh. From this figure, we can note that, as in the previous tests, approaches which employ the boundary DOFs (1)b are more robust with respect to the choice of the stabilization term and lead to more accurate results in terms of the pressure error (5.1) when the polynomial degree $k$ is not too high. Conversely, the Ortho (a) approach performs better than the Ortho (b) in terms of the velocity error (5.2).

The same conclusions can be reached by looking at Figures 29 and 30 that show the trend of errors as the mesh size h varies, for $k = 0$ and $k = 2$, respectively. Moreover, for the lowest polynomial degree $k = 0$, the locking phenomenon occurs and we can note it by observing a severe loss of accuracy mainly related to the velocity variable in the approaches Ortho (a)-S$_1$ and Ortho (a)-S$_4$. It should be noted that, in this test case, the right convergence rate is not restored even by increasing the polynomial degree $k$. For the higher values of the polynomial degree $k$ and for each tested approach, it is worth noting that both the errors (5.1) and (5.2) are mainly determined by local errors associated with very few cells near the point $(0, 0)$ where the diffusion tensor is not well-defined. These local errors are several orders of magnitude greater than the local errors characterizing the other cells. In particular, the pairs Ortho (a)-S$_1$ and Ortho (a)-S$_4$ exhibit very unstable behaviour near the point $(0, 0)$, as illustrated in Figure 31. This figure shows a comparison of the numerical solutions along the $y = x$ line for the approaches Ortho (a)-S$_1$, Ortho (b)-S$_1$ and Ortho (a)-S$_2$. We omit to show the numerical solution related to the approach Ortho (a)-S$_4$ since its behaviour is very similar to the one related to the approach Ortho (a)-S$_1$. This figure highlights how the choice of stabilization relies on the definition of the local degrees of freedom.

6.   Conclusions

In this paper, we carried out the analysis of the robustness of the mixed Virtual Element Method when problems characterized by highly anisotropic diffusion tensors are considered. Furthermore, a new set of boundary degrees of freedom based on moments computed against an $L^2([0, 1])$-orthonormal basis is also introduced.

Here, we report the results obtained on a set of benchmark problems by resorting to various approaches which differ for the sets of both the internal and the boundary degrees of freedom. For each benchmark problem, we propose different kinds of the stabilization term and we test the sensitivity of each proposed approach to the choice of the stabilization term in terms of both the condition number of the system matrix and of the errors (5.1) and (5.2).

In particular, the new set of boundary degrees of freedom seems to be more favourable in terms of errors by leading to a downward shift of the error curves, although, this choice generally does not ensure obtaining an improvement in the conditioning of ${\bf K}$. Indeed, the condition number of the system matrix seems to be mainly controlled by the choice of internal DOFs and by the anisotropic ratio.

Finally, the D-recipe version of the stabilization term with unit constant seems to be a good alternative to build a robust method for highly anisotropic diffusion problems.

Use of AI tools declaration

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

Acknowledgments

The author S.B. kindly acknowledges partial financial support provided by PRIN project "Advanced polyhedral discretisations of heterogeneous PDEs for multiphysics problems" (No. 20204LN5N5_003) and by PNRR M4C2 project of CN00000013 National Centre for HPC, Big Data and Quantum Computing (HPC) (CUP: E13C22000990001). The author S.S. kindly acknowledges partial financial support provided by INdAM-GNCS through project "Sviluppo ed analisi di Metodi agli Elementi Virtuali per processi accoppiati su geometrie complesse" and that this publication is part of the project NODES which has received funding from the MUR-M4C2 1.5 of PNRR with grant agreement no. ECS00000036. The author G.T. kindly acknowledges financial support provided by the MIUR programme "Programma Operativo Nazionale Ricerca e Innovazione 2014–2020" (CUP: E11B21006490005). Computational resources are partially supported by SmartData@polito. The authors are members of the Italian INdAM-GNCS research group.

Conflict of interest

The authors declare no conflicts of interest.