High-order interpolatory Serendipity Virtual Element Method for semilinear parabolic problems

We propose an efficient method for the numerical approximation of a general class of two dimensional semilinear parabolic problems on polygonal meshes. The proposed approach takes advantage of the properties of the serendipity version of the Virtual Element Method, which not only reduces the number of degrees of freedom compared to the original Virtual Element Method, but also allows the introduction of an approximation of the nonlinear term that is computable from the degrees of freedom of the discrete solution with a low computational cost, thus significantly improving the efficiency of the method. An error analysis for the semi-discrete formulation is carried out, and an optimal estimate for the error in the $L_2$-norm is obtained. The accuracy and efficiency of the proposed method when combined with a second order Strang operator splitting time discretization is illustrated in our numerical experiments, with approximations up to order $6$.


Introduction
In this work we present an interpolatory or quasi-interpolatory Serendipity Virtual Element Method (S-VEM) applied to semilinear parabolic equations on a space-time domain Q T = Ω × (0, T ), where Ω ⊂ R 2 is a polygonal domain and T > 0 The nonlinear function f : R → R is assumed to be globally Lipschitz continuous, i.e., there exists a constant L f > 0 such that the following bound holds The model (1.1) is found in many important applications such as: battery modeling [33], crystals growth [23], population dynamics [29], and in many other models in chemistry [34,27] and biology [25]. However, given the different nature of nonlinear terms, the task of finding exact solutions for such kind of problems becomes extremely demanding or even impossible. For that reason, there is a high interest in the development of efficient, accurate and robust numerical methods to approximate their solution. Since this work specifically concerns to the advantages of the serendipity version of the Virtual Element Method applied to the problem (1.1), an extensive list of numerical methods previously applied to this problem is out of our scope. The Virtual Element Method (VEM) is a novel technique for the numerical approximation of PDEs, introduced by Beirão da Veiga et. al., in [9], for an elliptic problem, and can be seen as a sensible extension of the classical finite element method to meshes with almost general polygonal elements. Discrete VEM spaces contain non-polynomial functions; however, such functions are not needed to be expressly known, as the discrete operators are computed through projections onto the space of piecewise polynomials of a given degree, which are computable using only some suitably chosen degrees of freedom (DoFs). Besides the advantages that come from the versatility of polygonal meshes, such as the natural use of "non-conforming" grids, more efficient and easier adaptivity and geometric approximation, robustness to mesh deformation, among others; the Virtual Element Method also allows for the imposition of conformity conditions on the global discrete spaces without struggling to explicitly compute their basis functions.
So far, the Virtual Element Method has been successfully applied to many important physical applications. In particular, recent efforts have been devoted to show the accuracy and advantages of this method in the numerical approximation of the solution to nonlinear problems such as: the Cahn-Hilliard equation [8,24], models in cardiology [7], nonlinear elasticity [18], the nonlinear Brinkman equation [22,28], bulksurface reaction-diffusion systems [21], pattern formation [19]; and semilinear elliptic [14,5], hyperbolic [2] and parabolic [1] equations.
In this work we aim to extend the idea of Adak and Natarajan in [3] to high-order approximations. In [3], the authors proposed a VEM discretization for the sine-Gordon equation with an interpolatory approximation of the nonlinear term, thus significantly reducing the computational cost of the method at each time step. Nevertheless, the main limitation of the technique in [3] is that it is only valid for approximations with k = 1, i.e., with the same order of convergence as polynomial approximations of degree k = 1. This is due to the fact that for k ≥ 2, the method requires some internal-moment DoFs, which unfortunately prevents a direct extension to high order approximations (see Remark 3.3 and Section 7 in [3] dedicated to discuss this limitation). Our idea to overcome this severe restriction relies on the serendipity variation of the VEM, introduced by Beirão da Veiga et. al., in [11], and later discussed by Russo in [30]. The main motivation of the S-VEM is indeed to reduce the number of internal-moment DoFs. Moreover, under certain conditions on the mesh, it is possible to completely eliminate them. It is also worth mentioning that the Serendipity VEM on quadrilaterals does not suffer from distortion as it is common for the serendipity FEM, see [11].
The main novelty and features of the proposed scheme are summarized as follows: a) To the best of our knowledge, this is the first time to use the S-VEM as spatial discretization for semilinear parabolic problems; for which the enhanced VEM has been preferred. b) An interpolant of the nonlinear term in the S-VEM space is introduced in the semi-discrete formulation. Under a condition on the degree of accuracy k that is associated to the geometric properties of the mesh, such interpolant is computed by simply evaluating the nonlinear term f (·) at the DoFs of the discrete solution. Thus significantly reducing the computational cost of the method, as it completely avoids the use of quadratures at each time step. c) An optimal error estimate for the semi-discrete formulation in the L 2 -norm is proven in spite of the use of such interpolant to approximate the nonlinear term. d) When the time variable is discretized by the symmetric Strang -operator splitting (SS-OS) time marching scheme, the nonlinear substeps can be decomposed as a set of completely independent one dimensional nonlinear problems, which makes the method naturally suitable for parallel implementations. e) In those elements of the mesh (if any) where the condition on the degree k is not satisfied, the interpolant of the nonlinear term is not computable from the DoFs of the discrete solution. In that case, we use a quasi-interpolatory approximation of the nonlinear term that also belongs to the local VEM space but is computable. Optimal error estimates and suitability for a parallel implementation are preserved.
The paper is structured as follows: in Section 2 we present the basic ideas and necessary projections for the description of the proposed method in the case when no internal-moment DoFs are needed, that throughout the paper we will refer to as "the ideal case". Optimal error estimates of order O(h k+1 ) in the L 2 -norm are proven for the S-VEM semi-discrete formulation in Section 3. In Section 4, an efficient fullydiscrete scheme, obtained by combining our interpolatory S-VEM discretization in space with a symmetric Strang -Operator Splitting time marching scheme is presented. The extension of the method to the general case when some internal-moment DoFs are needed, as well as the most important differences in the error estimate and the fully-discrete scheme are discussed in Section 5. Some numerical experiments, validating the accuracy and efficiency of the proposed method are included in Section 6. We end this work with some concluding remarks in Section 7.

Serendipity VEM discretization
Let T h be a polygonal partition of Ω and let h := max {h E | E ∈ T h } be the mesh size, where h E denotes the diameter of E. We first define, for each polygon E ∈ T h , the following enlarged local Virtual Element space [6]: where P k (e) and P k (E) denote the spaces of polynomials of degree at most k ≥ 1 on e and E, respectively. The DoFs uniquely identifying a function v ∈ V k (E) are choosen as the following linear functionals i) The values of v at the vertices of E.
ii) The values of v at the (k − 1) internal Gauss-Lobatto nodes on each edge e of E.
iii) The internal-moments: α=1 is a basis of P k (E) and r k := dim(P k (E)). The space V k (E) requires many more internal DoFs than the original VEM space presented in [9], but it readily provides enough information to compute the L 2 -projection of any v ∈ V k (E) onto P k (E). In practice, a subspace of V k (E) still containing all polynomials of degree at most k on E, is used as local VEM space. The basic idea to construct such subspace is to take the set of functions in V k (E) sharing some internal-moment DoFs with their projection onto the space P k (E); which gives origin to the so called, enhanced [6] and serendipity [11] versions of the VEM.
We first focus on the ideal case, where the Serendipity VEM does not require any internal-moment DoFs. An integer number η E ≥ 3 is associated to each element E ∈ T h , where η E is defined as the number of distinct straight lines containing at least one edge of E. In particular, if E is an N-sided strictly convex polygon without split edges, then η E = N. In the ideal case, the degree of accuracy k satisfies the condition: which in the spirit of the Serendipity VEM, allows the definition of a global discrete space whose associated DoFs are all node evaluations at the skeleton of the mesh, without requiring any internal-moment degree of freedom from the set iii).
The following projectors are needed to define the local S-VEM space and to present our semi-discrete formulation: • The Ritz-Galerkin projection π ∇ k,E : Using the Green's formula, the projection π ∇ k,E (·) is computable for any v ∈ V k (E) using the degrees of freedom i), ii) and iii), see [9,Sect. 4.5].
• The standard L 2 -orthogonal projector π 0 k,E : which is directly computable from the set of DoFs iii).
• The "boundary" projector π ∂ k,E : that is well-defined under the condition k < η E , and can be computed using only the boundary DoFs in sets i) and ii).
In the ideal case, the local S-VEM space that only requires the boundary DoFs from sets i) and ii). Therefore, for an N-sided polygon E, dim V S k (E) = Nk. The global S-VEM space is consequently defined as A representation of the DoFs for the local space of the original VEM in [9], for different N-sided polygons and the maximum degree k satisfying the aforementioned condition is presented in Fig. 1. We emphasize that, in all these cases, the local space V S k (E) does not require any of the internal-moment DoFs represented by red squares in the figure.
As mentioned before, any function v ∈ V S k (E) is uniquely determined by its boundary DoFs. Denoting by d S k,E = dim(V S k (E)), and numbering the nodes associated to the DoFs as ξ i , with i = 1, . . . , d S k,E , we can define the linear functionals dof i : A natural basis arise, by taking the canonical basis functions The following interpolatory representation is then obtained for such representation, in turn, allows us to define the interpolant operator I k h : C 0 Ω → V S k (T h ) whose restriction to each element E ∈ T h is defined as follows: for all g ∈ C 0 E ,

Semi-discrete formulation
The weak formulation of the model problem (1.1) is: find u ∈ L 2 0, T, Analogously, our semi-discrete interpolatory S-VEM formulation seeks an approximation where the bilinear forms a h : For each element E ∈ T h , the restrictions a E h (·; ·) and m E h (·; ·) are split into a consistency and a stability parts by where I denotes the identity operator and the stabilization terms s E a (·; ·), s E m (·; ·) are symmetric bilinear forms scaling as a E (·; ·) and m E (·; ·), respectively; more precisely, there exist positive constants α 1 , α 2 , β 1 , β 2 such that In fact, there are many possible choices for the stability terms; however, in our implementation we will limit ourselves to use a very simple stabilization proposed in [9], namely, the dofi-dofi. For a thorough study of different stability choices, see [12,26]. By construction, both a h (·; ·) and m h (·; ·) satisfy the following two important conditions: • Stability: There exist mesh-independent positive constants α * , α * , β * , β * such that The last term in the semi-discrete variational formulation (2.6a) satisfies the following crucial identity which clearly shows that, in the ideal case, the computation of the nonlinear term requires only the matrix representation of the bilinear form m h (·; ·) and the evaluation of f (·) at the DoFs of the discrete solution u h .
Remark 1. Applying stabilization in the last term of (2.6a) in the semi-discrete formulation is not necessary to obtain optimal convergence, but it can be computationally convenient, as shown in Section 4.

Remark 2. The initial condition approximation u
is suitable for imposing random initial data, which is commonly of interest in this kind of problems.

Error analysis
This section is devoted to get an optimal error estimate in the L 2 -norm for the solution of the semidiscrete formulation (2.6). The main ideas are taken from the error analysis carried out in [35] for the enhanced VEM applied to linear parabolic problems and its recent extensions to semilinear parabolic problems [5,1,4]. Nonetheless, in Theorem 1 we address the following differences: • The approximated solution is sought in the S-VEM space V S k (T h ).
• The nonlinear term is approximated by its ; v h in the semi-discrete formulation (2.6) includes also a stabilization part, that was not present in the formulation in [3] for the sine-Gordon equation. This choice is computationally convenient when using an operator splitting time marching scheme, while retaining the same optimal convergence.
• Pure homogeneous Neumann boundary conditions are considered.
In what follows we will make the following assumptions on the mesh: There exists a constant ρ > 0, such that every element E ∈ T h is star-shaped with respect to a ball B := B ρh E (z) centered at z ∈ E and with radius ρh E , where h E := diam (E). In addition, every edge e of E satisfies |e| ≥ ρh E .
The above assumption guarantees that the following condition holds: for each E ∈ T h , there exists a "virtual triangulation" T E of E such that T E is uniformly shape regular and quasi-uniform. The corresponding mesh size of the auxiliary triangulation T E is proportional to h E = diam (E) and each edge of E is a side of a triangle in T E .
The elliptic projection operator P h : Since P h (u) is the solution to the variational problem (3.1), by the coercivity and continuity of a h (·; ·) and the continuity of the linear functional a (u; ·), the projection operator P h is well-defined. Furthermore, we can prove the following estimate as in [35, Lemma 3.1].

Lemma 1.
Let Ω be a convex domain, and u ∈ H k+1 (Ω). Under Assumption 1, there exists a constant C α > 0, depending on α * and α * in (2.9a) but independent of h such that Using standard arguments as in [9] and the classical Dupont-Scott theory in [13], the following estimates for the interpolant I k h (·) and the projection π 0 k (·) are obtained.
, there exists a polynomial u π ∈ P k (E), and a positive constant C π , depending only on k and ρ, such that In Lemma 4, a norm equivalence for the S-VEM space is introduced in order to exploit the Lipschitz property of f (·) in the analysis.
. Lemma 4 is an extension of the classical results in [20] for finite element spaces and can be proven following the arguments used by Chen and Huang in [17, Thm. 4.5 and Corollary 4.6].
From the above Lemma we can derive the following important bound in our error analysis.
Lemma 5. The following bound holds for any u ∈ C 0 E and u h ∈ V S k (E)
Proof. By Lemma 4, and the Lipschitz continuity of f (·) we have To conclude the proof it suffices to use the triangle inequality in the last term.
The following theorem provides the optimal error estimate for the semi-discrete formulation (2.6) under suitable regularity conditions for the exact solution. We will use C to denote a generic constant independent of the mesh size h and the arguments of the functions in the proof will be omitted unless they are necessary. (3.7) Proof. We start decomposing e u := u − u h as e u = ξ u − θ h , where ξ u = u − P h (u) and θ h = u h − P h (u). From Lemma 1 and the identity u(·, t) = u(·, 0) + t 0 u t (·, τ)dτ we have the following bound for ξ u Therefore, in order to get the desired estimate, it only remains to bound θ h (·, t) L 2 (Ω) . We now proceed similarly as in [1]. Since θ h ∈ V S k (T h ), adding and subtracting appropriate terms in the semi-discrete formulation (2.6a), for any hence, we will look for local estimates of T E 1 and T E 2 on each element E ∈ T h . By the k-polynomial consistency property (2.8b) of m h (·; ·), we can decompose T E 1 as (3.10) By the Cauchy-Schwarz inequality and Lemma 3, it is easy to see that On the other hand, by the stability of the L 2 -orthogonal projection π 0 k (·) and Lemma 2 we have To bound R 3 , we first observe that by the triangle inequality, the stability of the L 2 -orthogonal projection and Lemmas 2 and 3 we have (3.13) Lemma 5 and the bound (3.13) together with the triangle inequality and the continuity of m h (·; ·) provide the following estimate for R 3 : In a similar way, decomposing and applying similar steps as before, by the commutativity of ∂ ∂t (·) and P h (·), we get the following bound Integrating from 0 to t at both sides of (3.9) and taking v h = θ h ; since a h (θ h ; θ h ) ≥ 0, by the estimate (3.8) for ξ u , the bounds (3.11)-(3.15) and Young's inequality, we get the following estimate + u t 2 L 1( 0,t,H k+1 (Ω)) + u t 2 L 2( 0,t,H k+1 (Ω)) + u 2 L 2( 0,t,H k+1 (Ω)) + f (u) 2 L 2( 0,t,H k+1 (Ω)) , and by Grönwall's lemma, since θ h (·, 0) = P h (u 0 ) − u 0 + u 0 − I k h u 0 , combined with the bound (3.8) for ξ u , and the estimates (3.2)-(3.3), we get the desired estimate (3.7) in our theorem.
in ( On the other hand, bound (3.13) is not necessary when the stability part of the last term of (2.6a) in the semi-discrete formulation is not considered.

Fully-discrete scheme
It is evident that the efficiency of any ODE solver applied to (2.6) will be greatly benefited from the fast evaluation of the nonlinear term in our semi-discrete formulation. In this paper, we choose the second order symmetric Strang operator splitting (SS-OS) method [32] as time marching scheme to illustrate the advantages of the proposed technique.
Denoting by M and A the matrix representation of the bilinear forms m h (·; ·) and a h (·; ·), respectively; by the identity (2.10), the semi-discrete formulation (2.6) can be written as a system of nonlinear differential equations as where U h is the vector of the representation coefficients of u h in the basis of V S k (T h ); and the components of the vector In the ideal case, f h (U h ) is the vector obtained from a component-wise evaluation of the nonlinear function f (·) at the entries of U h .
The SS-OS time marching scheme decomposes the system of differential equations (4.1) as a series of linear and nonlinear substeps, usually associated to diffusion and reaction terms, of the form DRD decomposition: RDR decomposition: where τ = t n+1 − t n and U n h is the vector approximation of u h (·, t n ). The efficiency of combining some discontinuous Galerkin methods with an interpolatory approximation of the nonlinear term as spatial discretization on classical meshes with the SS-OS time marching scheme was assessed by Castillo and Gómez in [15,16].
A necessary condition to retain the second order accuracy of the full SS-OS step is that each substep in (4.2) must be solved with a second order ODE solver itself. Although we are free to choose the solver for each step, implicit methods might be more appropriate. Conversely, if an explicit method were used, we would face a very restrictive CFL condition associated to the linear substeps, while for the nonlinear substeps the method might become unstable in the case of stiff nonlinearities.
From the discussion above we decide to apply the Crank-Nicolson method to each substep in (4.2). For the DRD decomposition (4.2a) the resulting fully-discrete method reads The following remarks are in order: • The linear substeps (4.3a) and (4.3c) only consist in solving two linear systems with the same matrix. For a fixed time step τ such matrix is even the same at any time, which is advantageous since a preconditioner or a full Cholesky factorization can be computed just once at the beginning of the simulation.
• The nonlinear substep requires the solution of the nonlinear system (4.3b), which is completely independent for each component of the vector U (2) h , and as such, highly parallelizable. Note that we have cancelled matrix M at both sides of this equation. Such cancellation is only possible because stabilization was also applied to the nonlinear term in (2.6a); otherwise, a large coupled system of nonlinear equations would be obtained. If we apply the Newton's method to (4.3b) each linear iteration s reads is also diagonal, the solution of (4.4a) reduces to a trivial entry-by-entry division.
We end this section with the following well-posedness result of the fully-discrete scheme. Proposition 1. The fully-discrete schemes DRD and RDR are well-posed for any 0 < τ < 2/L f .
Proof. Without loss of generality we will prove the well-posedness only for the DRD scheme. Since matrix M + τ 4 A is symmetric and positive definite, the existence of the solution of each linear substep in (4.3a) and (4.3c) is guaranteed.
On the other hand, each independent one dimensional problem in the nonlinear substeps (4.3b) is equivalent to find a fixed point of the function g(x) = a − τ 2 f (x) for some constant a, which can be easily shown to be a contraction as long as 0 < τ < 2/L f ; therefore, under such condition, the existence of the solution of the nonlinear substeps is also guaranteed.
Existence and uniqueness of the full step in (4.3) would then proceed from those of each susbtep.

Extension to arbitrary k
We now present an extension of the interpolatory S-VEM to the general case, when some internalmoment DoFs are needed. The main drawback in such case is that for k ≥ η E , condition (2.1) is not enough to define a projection due to the existence of P k (E)-bubbles. Hence, some additional internal-moment DoFs and a computable projection operator are needed.
Since for non-convex polygons the choice of the additional DoFs is more involved, see [11,Sect. 3], we will focus on the case of convex polygons. For convex polygons, if the internal-moment DoFs up to order k − η E are added, the projection π S k,E : V k (E) → P k (E) defined in [30] as is well-defined and computable from the DoFs by definition. For any convex polygon E ∈ T h , if k ≥ η E , the local Serendipity VEM space is then defined as

Numerical experiments
In this section we present some numerical experiments to show the accuracy and efficiency of the proposed scheme. An object oriented implementation in MATLAB was developed for high order approximations on general polygonal meshes. As time marching scheme we use the SS-OS method (4.2) presented in section 4. All the linear systems were solved with the preconditioned conjugate gradient (PCG) method. The incomplete Cholesky factorization with a drop tolerance of 10 −5 was used as preconditioner. Linear and nonlinear systems were solved with a tolerance of 10 −10 as stopping criteria; and numerical quadratures for each polygon were obtained using the Vianello approach [31]. The sets of meshes used in all the experiments are exemplified in Fig. 2. Note that, for these meshes, the values of η E can be known a-priori, as for strictly convex N-sided polygons η E = N and all the non-convex polygons in Fig. 2c satisfy η E ≥ 8.  In order to illustrate the accuracy and efficiency of the proposed method, we will compare our results with those obtained for the enhanced VEM proposed in [1]. While the linear substeps of the SS-OS time marching scheme are similar for both versions, the nonlinear substeps for the method in [1] require to solve the following strongly coupled system of nonlinear equations where F h (·) is the nonlinear operator defined as The nonlinear systems (6.1) will be solved using a semilinear iterative method, that avoids computing the Jacobian of the nonlinear term. Each linear iteration s consists in solving the following linear system On the other hand, the reported execution times correspond to computations carried out on a DELL laptop with an Intel Core i7-8750h processor, 32Gb of RAM and Linux operating system.

Accuracy test
As first experiment we numerically asses the accuracy of the proposed method. We consider a manufactured problem on Q T = (0, 1) 2 × (0, 1] with a nonlinear term f (u) = 1/(1 + u 2 ), adding a source term such that the exact solution be u(x, y, t) = e −t cos(πx) cos(πy).
In Fig. 3 we present the errors in the L 2 -norm with respect to π 0 k (u h ) at the final time T , i.e., at Σ T := Ω× {T }, for each kind of mesh. In the same plot we have included the errors obtained by the enhanced VEM in [1] as reference; and no significant difference in terms of accuracy is observed. Optimal rates of convergence of order O h k+1 are obtained as stated in Theorem 1. The time step was taken as τ = O h (k+1)/2 in order to equilibrate the errors in space and time.  (c) Non-convex meshes. To evaluate the temporal accuracy of the fully-discrete scheme, we use a sequence of time refinements with τ = 1.25 × 10 −1 , 6.25 × 10 −2 , 3.125 × 10 −2 , 1.5625 × 10 −2 ; and in order to let the time error dominate, computations were carried out for the finest voronoi mesh and k = 4. The obtained rates of convergence for the DRD and the RDR splitting methods are shown in Fig. 4 and validate the second order in time O τ 2 accuracy of the SS-OS fully-discrete scheme. In this experiment, better accuracy is observed for the RDR splitting. Not shown here, similar results were obtained for the other meshes.
In Table 1, we compare the number of global degrees of freedom for the S-VEM and the enhanced VEM in [1], where naturally the reduction in the number of degrees of freedom depends on the mesh and a more noticeable reduction is obtained at increasing k. This is also illustrated in Fig. 5, where we compare the accuracy of both methods with respect to the number of DoFs.

Efficiency test
In this experiment we consider the following Allen-Cahn equation on Q T = (0, 1) 2 × (0, 22.5] as in [1]: u(x, y, 0) = cos(2πx 2 ) cos(2πy 2 ), in Ω, (a) DRD.    where the nonlinear term f (u) = u 3 − u only satisfies a local Lipschitz condition. In fact, the error estimate in Theorem 1 is still valid if f (·) is only locally Lipschitz continuous under the additional assumption of both the exact and the approximated solutions to be bounded. In order to show the efficiency of the proposed method, we compare our results with those obtained for the interpolatory VEM in [3] and the enhanced VEM in [1]. In all these experiments, we consider the finest meshes of each kind, τ = 5 × 10 −3 as time step and the RDR splitting.
In Table 2 we report the CPU execution times for the approximation of the Allen-Cahn equation (6.2) with ǫ = 0.01 for the proposed method and the interpolatory VEM presented in [3]. We recall that the method in [3] is limited to k = 1 and does not include the stability part of the nonlinear term, so the nonlinear systems in the SS-OS fully-discrete scheme (4.2) remain coupled. We observe that the times in the linear substeps are approximately equal in both cases, which is expected as both methods have the same number of DoFs. However, for the nonlinear substeps our method performs about 20 to 70 times faster depending on the mesh; and a total boost of approximately 10 times is obtained in all the cases.  2) in the test problem 6.2 with ǫ = 0.01, for the proposed S-VEM and the interpolatory VEM in [3] with k = 1.
In a similar way, in Table 3 we compare the CPU execution times for the proposed method and the enhanced VEM in [1] with different degrees of approximation. Since the proposed method requires less DoFs, it performs faster for the linear substeps. As for the nonlinear substeps, our method performs from 40 to 2500 times faster depending on the mesh and the degree of accuracy. A total boost of about 12 to 110 times is obtained. For each mesh, we have indicated those degrees where some internal moment DoFs are needed; in such cases, the extended version from Section 5 was used and a significant improvement in the efficiency of the method is still observed. The substantial reduction obtained for the non-convex mesh is a consequence of the high number of quadrature points required for the VEM in [1] Table 3: CPU execution times for the Allen-Cahn equation (6.2) in the test problem 6.2 with ǫ = 0.01, for the proposed S-VEM and the enhanced VEM in [1]. The red star symbol ⋆ indicates that the extended version from Section 5 was needed.
In Figure 6 we show the evolution of the approximated solution π 0 k (u h ) for the Allen-Cahn equation with ǫ = 0.01, which is expected to converge to its stable state u = −1. The plots portray the same behaviour observed in [1] for the enhanced VEM.

Conclusions
In this work, an interpolatory Serendipity Virtual Element method for semilinear parabolic problems on polygonal meshes is proposed. A significant reduction in the computational cost of the method is obtained by approximating the nonlinear term with an element of the S-VEM space. Optimal error estimates of order O h k+1 in the L 2 -norm are proven for the semi-discrete formulation.
To exploit the structure of the system of nonlinear differential equations arising from the semi-discrete formulation, we use a second order operator splitting as time marching scheme, which decouples the linear and nonlinear terms. In the ideal case, with only boundary DoFs, the nonlinear substeps consist in solving a set of completely independent one dimensional nonlinear equations; while in the extension proposed to the case when some internal-moment DoFs are required, it is also necessary to solve an additional set of independent small nonlinear systems on each element of the mesh that does not satisfy the condition of the ideal case. Our numerical experiments validate the optimal convergence of the method and the improvement in efficiency respect to the enhanced VEM in [1].