Virtual Element Approximation of Two-Dimensional
Parabolic Variational Inequalities

Variational inequalities have been an active research field in the last decades and has found many important applications in finance and engineering [50, 10, 57]. For example, they are used in the formulation of the one-phase Stefan problem [58]. The Allen—Cahn equation, one of the models of the kinetics of grain growth in polycrystals, can be treated as a parabolic variational inequality [23]. The American put option problem [49] becomes a one-phase Stefan problem after a suitable change of variable [48]. The electrochemical machine problem is also modeled using variational inequalities [39]. Static contact problems, frictional contact problems, and thermal expansion problems can be described using variational inequalities, cf. [32]. The numerical approximation of the solution to variational inequalities has also been a challenging area of research since both the design of numerical methods and the convergence analysis are not straightforward [45].

The Galerkin approach and, in particular, the finite element method (FEM) has proven to be quite effective to this purpose. The linear Galerkin FEM for the time-dependent parabolic variational inequality (with zero obstacle) was originally proposed in [47]. In this paper, which is the most pertinent to our current work, a priori error estimates in the $L^{\infty}$ norm are derived assuming that the solution is in $L^{\infty}\big{(}0,T;W^{2,p}(\Omega)\big{)}$ and its first derivative in time is in $L^{\infty}\big{(}0,T;H^{1}_{0}(\Omega)\big{)}\cap L^{\infty}\big{(}0,T;L^{\infty}(\Omega)\big{)}$ (we explain this notation and provide a formal definition of these functional spaces later in this section). A priori estimates in the $L^{2}$ norm are also derived for the Galerkin method in [41] assuming that the solution is in $L^{2}(0,T;L^{2}(\Omega))$ and under certain regularity assumptions on the angles of each element of the triangulations. In [64], $L^{2}$ error estimates are derived for a fully discrete scheme based on the $\theta$-method in time. Inspired by the American put option problem, a posteriori error estimates are studied in [53]. In [19], the Authors derive error estimates for the parabolic variational inequality problem in the uniform norm. Moreover, in [21, 46, 3, 59], mathematical models of the parabolic obstacle problem related to the American put option problem and the Stefan problem are investigated.

In this work, we consider the approach that was originally proposed in [47] for solving a parabolic variational inequality on triangular meshes and study how to generalize it to polygonal meshes using the virtual element method (VEM) [11, 2]. Designing Galerkin schemes for meshes with polygonal elements in 2D and polyhedral elements in 3D has been a major topic in the numerical literature of partial differential equations of the last two decades. Several classes of numerical methods have been designed that are suitable to meshes with elements having very general geometric shapes. Other than the VEM, a surely nonexhaustive list includes the polygonal/polyhedral finite element method (PFEM) [60, 61, 22], the mimetic finite difference method (MFD) [51, 14, 13], the hybridizable discontinuous Galerkin (HDG) method and the hybrid high-order (HHO) method [34, 38, 36, 37]. Pertinent to the topic of our work are also the papers of References [7, 6].

The virtual element method was proposed as a variational reformulation of the mimetic finite difference method of References [26, 13] for the Poisson equation, and later extended to the numerical approximation of general elliptic equations [12], elasticity problems [54], eigenvalue problems [55, 56, 33, 42], Stokes and Navier-Stokes equations [4, 28, 43, 15, 18], and the Cahn-Hilliard equations [5]. Furthermore, the mixed formulation [27] the nonconforming formulation [9, 30, 31], and the enriched formulation [20] have been proposed and a posteriori error estimations [16, 29, 17, 35] have been derived for mesh adaptivity. VEM for anisotropic polygonal discretizations are also found in [8].

The VEM satisfies a Galerkin-type orthogonalization property on polynomial subspaces and can be seen as a generalization of the FEM on arbitrary polytopal meshes. The finite dimensional approximation spaces consist of polynomial and nonpolynomial functions that satisfy a partial differential equation locally defined on the mesh elements. The nonpolynomial functions are not known inside the elements, but the degrees of freedom of the virtual element functions are carefully chosen so that some polynomial projection operators are computable. These projection operators make it possible to design computable bilinear forms for the discrete variational formulation. Since an explicit knowledge of the virtual element functions is not required in the practical implementation, such “virtual” setting works for very general shaped polytopal elements. For example, nonconvex elements and elements with hanging nodes are admissible and the latter do not require any special treatment.

Due to the expected low regularity of the solution, our method is based on the lowest-order approximation space proposed in [11, 2]. The degrees of freedom are the vertex values and our VEM coincides with the FEM of Reference [47] on all triangular meshes. The generalization to the virtual element framework is nontrivial and the design of an effective VEM and its analysis is challenging for several reasons that we illustrate below. First, the variational formulation is given on the subset of the nonnegative virtual element functions, which we identify with those functions of the virtual element space whose degrees of freedom, i.e., the vertex values, are nonnegative. The property that a function with nonnegative vertex values is nonnegative is obvious for a linear polynomial interpolating such values on a triangular element. However, to prove that such property holds for a virtual element function on a polygonal element is a nontrivial task. In fact, such functions are not generally known in closed form, but only as the solutions of an elliptic partial differential equation that is locally set on the polygonal element. We address this issue by noting that the lowest-order virtual element space consists of functions that are harmonic inside each element and have a continuous piecewise linear trace on the elemental boundary given by the interpolation of the vertex values. Consequently, we can prove the nonnegativity property by invoking the Maximum and Minimum Principle Theorem [44]. According to this theorem, a nonconstant harmonic function on a compact set of points, e.g., a (closed) polygonal element, must take its maximum and minimum value on the boundary. If all its vertex values are nonnegative, so is their piecewise linear interpolation on the elemental boundary and the function itself inside the element. Unfortunately, we cannot apply this theorem to the modified (“enhanced”) virtual element space introduced in [2] as its functions are no longer harmonic. This fact poses a major issue to the design of our VEM since we need an $L^{2}$-like orthogonal projector for the calculation of the mass matrix in the discretization of the time derivative term. To address this issue, we design a different projector, which is still computable from the degrees of freedom of the space and is orthogonal with respect to an approximate $L^{2}$ inner product. We carefully characterize the approximation properties of this operator to prove the convergence of the VEM, estimate the approximation error and derive the convergence rate for the refinement in time and space.

We also prove the well-posedness of the numerical method, i.e., existence and uniqueness of the virtual element solution, in two different ways. The first proof reveals the contractive nature of the scheme, which motivates an iterative implementation at every time step from a practical viewpoint. The second proof generalizes a minimization argument briefly mentioned in [47] to the new virtual element framework proposed in this work and establishes a clear connection between the VEM and the minimization of quadratic functionals.

To carry out the theoretical analysis and prove the convergence of the VEM, we investigate how the virtual element reformulation impacts on the original convergence proof of Reference [47]. A major ingredient of the latter is the existence of a nonnegative quasi-interpolation operator for functions that are only $H^{1}$-regular as, for example, the derivative in time of the parabolic inequality solution. To address this point, we generalize the construction of such operator in [47], so that it can work on polygonal elements with the desired nonnegativity property. Finally, we identify the new terms that arise from the “variational crime” introduced by the virtual element method and provide an upper bound for all of them.

The numerical experiments confirm the validity of our approach by solving a manufactured solution problem on very general meshes including distorted square elements, nonconvex elements and Voronoi tesselations. The experimental convergence rates reflects the convergence rates expected from the theoretical analysis.

The outline of the paper is as follows. In the rest of this section, we introduce some background material from functional analysis and the notation used in the paper. In Section 2, we discuss the continuous weak formulation of the mathematical model. In Section 3, we present our virtual element method for the parabolic inequality problem. In Section 4, we introduce some technical lemmas and detail the construction of the quasi-interpolation nonnegative operator for the convergence analysis. In Section 5, we prove the convergence of the method and derive the a priori error estimate. In Section 6, we assess the performance of the method on three different families of polygonal meshes. In Section 7, we summarize our results and offer the final remarks.

We let $\mathcal{K}=\big{\{}\,v\in H^{1}_{0}(\Omega):\,v\geq 0~{}\textrm{a.e.~{}in~{}}\Omega\,\big{\}}$ be the subset of the nonnegative functions in $H^{1}_{0}(\Omega)$. We also consider the positive real number $T$ representing the final integration time and the time interval $J=[0,T]$, and introduce the bilinear form

This bilinear form is coercive and continuous on $H^{1}_{0}(\Omega)$. So, there exists two real, positive constants $\alpha$ and $M$ such that that $\alpha\left\|v\right\|_{1}^{2}\leq a(v,v)$ and $a(v,w)\leq M\left\|v\right\|_{1}\left\|w\right\|_{1}$ for all $v$, $w$ in $H^{1}_{0}(\Omega)$. We search the solution $u(t)$ to the parabolic variational inequality problem for a given right-hand side source term $f$ and initial state $u_{0}$, which reads as

Find $u(t):J\to\mathcal{K}$ such that, for almost every $t\in J$ it holds that

The solution $u$ exists and is unique [25] under the assumptions

(A1) $f\in C{}\big{(}J;L^{\infty}(\Omega)\big{)}$;

(A2) $\partial f/{\partial t}\in L^{2}\big{(}J;L^{\infty}(\Omega)\big{)}$;

(A3) $u_{0}\in W^{2,\infty}(\Omega)\cap\mathcal{K}$.

In particular, if assumptions (A1)-(A3) are true, solution $u$ is such that:

where $\partial^{+}u/{\partial t}$ denotes the right-hand derivative of $u$ with respect to $t$. Moreover, $u$ satisfies the partial differential equations

where, for almost every $t\in J$, $\Omega^{+}(t)=\big{\{}\mathbf{x}\in\Omega:u(\mathbf{x},t)>0\big{\}}$ and $\Omega^{0}(t)=\big{\{}\mathbf{x}\in\Omega:u(\mathbf{x},t)=0\big{\}}$.

Finally, we partition the time interval $[0,T]$ into $N$ equally spaced subintervals $J_{n}=\big{[}t^{n},t^{n+1}\big{]}$ having size $\Delta t=t^{n+1}-t^{n}=T/{N}$, and let $m(\Gamma_{n})$ denote the area of the set

Our last assumption is that

(A4) $~{}\sum_{n=0}^{N-1}m(\Gamma_{n})\leq C$ for some real, positive constant $C$ independent of $h$ and $\Delta t$.

This assumption together with (A1)-(A3) will be used in the convergence analysis of the method that we perform in Section 5.

Let $V_{h}$ be a conforming finite dimensional subspace of $H^{1}_{0}(\Omega)$ that will be referred to as the virtual element space. Let $m_{h}(\cdot,\cdot),\,a_{h}(\cdot,\cdot):V_{h}\times V_{h}\to\mathbbm{R}$ be the virtual element approximation of the $L^{2}$ inner product $(\cdot,\cdot)$ and the bilinear form $a(\cdot,\cdot)$. Let $f_{h}$ be the element of $(V_{h})^{\prime}$, the dual space of $V_{h}$, such that $(f_{h},\cdot):V_{h}\to\mathbbm{R}$ is a virtual element approximation of the linear functional $(f,\cdot)$ (we use the same symbol $f_{h}$ to denote the Ritz representative of $f_{h}$ in $V_{h}$). Then, we introduce the finite-dimensional subset $\mathcal{K}_{h}=V_{h}\cap\mathcal{K}=\big{\{}v_{h}\in V_{h}:v_{h}\geq 0\textrm{~{}in~{}}\Omega\big{\}}$ of the virtual element functions that are nonnegative in $\Omega$. We denote the evaluation of a time-dependent quantity $w(t)$ at $t^{n}$ by $w^{n}=w(t^{n})$, and define the discrete difference operator $\partial w^{n}=\big{(}w^{n+1}-w^{n}\big{)}/{\Delta t}$, which provides the time variation of $\{w(t^{n})\}_{n}$ in the time interval $\big{[}t^{n},t^{n+1}\big{]}$.

The virtual element approximation $U_{h}^{n}$ to $u(t^{n})$ is the solution of the following discrete problem: Find $\{U_{h}^{n}\}_{n=0,\ldots,N}$ with $U_{h}^{n}\in\mathcal{K}_{h}$ for every $n=0,1,\ldots,N$ such that

for every $v_{h}\in\mathcal{K}_{h}$ with the initial solution field $U_{h}^{0}$ satisfying

This section is devoted to the definition of $V_{h}$, the construction of the bilinear forms $m_{h}(\cdot,\cdot)$ and $a_{h}(\cdot,\cdot)$ and the linear functional $(f_{h},\cdot)$, and the characterization of their approximation properties. Furthermore, a possible choice of the initial approximation of $u_{0}$, i.e., $U_{h}^{0}$, which satisfies (13), is provided by could be chosen as $U_{h}^{0}=I_{h}u_{0}$, where $I_{h}$ is the interpolation operator that will be defined in Section 4.2.