Auxiliary space preconditioners for virtual element methods on polytopal meshes

TL;DR

This paper introduces auxiliary space preconditioners for virtual element methods on polytopal meshes, improving the efficiency of solving linear systems from second order elliptic equations with robust performance regardless of mesh size or coefficient jumps.

Contribution

The paper develops novel auxiliary space preconditioners specifically designed for virtual element methods on polytopal meshes, ensuring uniform boundedness of condition numbers.

Findings

Condition numbers are uniformly bounded, independent of problem size.

Preconditioners perform well across various numerical experiments.

Effective for second order elliptic equations with coefficient jumps.

Abstract

In this paper, we develop the auxiliary space preconditioners for solving the linear system arising from the virtual element methods discretization on polytopal meshes for the second order elliptic equations. The preconditioners are constructed based on an auxiliary simplicial mesh. The condition numbers of the preconditioned systems are uniformly bounded, independent of the problem size and the jump in coefficients. Several numerical experiments are presented to demonstrate the performance of the preconditioners.

Tables

Table 1. Table 1. Estimated condition numbers (number of PCG iterations) in 2D with constant coefficients.

# Polytopal Elements	10	${10}^2$	${10}^3$	${10}^4$	${10}^5$
$𝒦\left(A\right)$	3.45 (9)	3.86e1 (41)	3.80e2 (117)	3.88e3 (351)	4.07e4 (1100)
$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.07(6)	3.78 (15)	3.20e1 (37)	3.17e2 (104)	3.17e3 (318)
$𝒦\left(B_{\mathrm{fict}}A\right)$	2.92 (8)	5.75 (26)	7.53 (29)	8.73 (32)	9.67(36)
$𝒦\left(B_{\mathrm{add}}A\right)$	1.53 (9)	1.71 (14)	1.94 (14)	1.99 (14)	2.00 (13)
$𝒦\left(B_{\mathrm{mul}}A\right)$	1.06 (8)	1.21 (10)	1.04 (7)	1.02 (6)	1.02 (6)

Table 2. Table 2. Estimated condition numbers (number of PCG iterations) in 2D with jump coefficients.

# Polytopal Elements	10	${10}^2$	${10}^3$	${10}^4$	${10}^5$
$𝒦\left(A\right)$	2.44 (11)	2.73e6 (578)	-	-	-
$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.18(5)	3.90e2 (26)	3.93e3 (409)	-	-
$𝒦\left(B_{\mathrm{fict}}A\right)$	3.27 (8)	6.94 (33)	6.42 (36)	11.6 (44)	13.6 (53)
$𝒦\left(B_{\mathrm{add}}A\right)$	1.54 (9)	3.51 (20)	3.60 (25)	3.67 (25)	3.80 (26)
$𝒦\left(B_{\mathrm{mul}}A\right)$	1.06 (6)	1.74 (15)	1.82 (16)	1.84 (16)	1.88 (17)

Table 3. Table 3. Estimated condition numbers (number of PCG iterations) in 2D voronoi polygonal mesh.

# Polytopal Elements	10	${10}^2$	${10}^3$	${10}^4$	${10}^5$
$𝒦\left(A\right)$	4.76 (9)	6.89e1 (52)	6.59e2 (171)	6.49e3 (537)	-
$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.13 (6)	4.93 (17)	3.81e1 (45)	3.57e2 (134)	3.40e3 (400)
$𝒦\left(B_{\mathrm{fict}}A\right)$	4.66 (9)	7.92 (34)	2.04e1 (43)	2.32e1 (46)	1.62e1 (52)
$𝒦\left(B_{\mathrm{add}}A\right)$	1.58 (9)	1.72 (16)	3.09 (18)	3.16 (19)	1.91 (17)
$𝒦\left(B_{\mathrm{mul}}A\right)$	1.32 (11)	2.25 (16)	1.48 (13)	1.29 (12)	1.14 (10)

Table 4. Table 4. Estimated condition numbers (number of PCG iterations) in 3D. The coefficient κ | Ω 1 ∪ Ω 2 = κ 1 = 10 k evaluated-at 𝜅 subscript Ω 1 subscript Ω 2 subscript 𝜅 1 superscript 10 𝑘 \kappa|_{\Omega_{1}\cup\Omega_{2}}=\kappa_{1}=10^{k} for various choices of k 𝑘 k , and κ | Ω ∖ ( Ω 1 ∪ Ω 2 ) = 1 . evaluated-at 𝜅 Ω subscript Ω 1 subscript Ω 2 1 \kappa|_{\Omega\setminus(\Omega_{1}\cup\Omega_{2})}=1.

${\kappa }_1$	$h$	$2^{-2}$	$2^{-3}$	$2^{-4}$	$2^{-5}$	$2^{-6}$	$2^{-7}$
${10}^{-6}$	$𝒦\left(A\right)$	1.15e6 (8)	8.76e6 (28)	6.94e7 (56)	5.54e8 (110)	4.43e9 (215)	3.54e10 (420)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (23)	2.44e1 (43)	9.57e1 (71)	3.81e2 (118)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.44 (11)	1.41 (10)	1.39 (9)	1.37 (8)	1.33 (7)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (10)	1.35 (11)	1.73 (15)	1.92 (17)	1.98 (17)	1.99 (16)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)
${10}^{-4}$	$𝒦\left(A\right)$	1.15e4 (7)	8.76e4 (26)	6.94e5 (51)	5.54e6 (99)	4.43e7 (194)	3.54e8 (379)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (23)	2.44e1 (38)	9.57e1 (63)	3.81e2 (115)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.44 (11)	1.41 (10)	1.39 (9)	1.37 (8)	1.33 (7)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (9)	1.35 (11)	1.73 (14)	1.92 (15)	1.98 (15)	1.99 (13)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)
${10}^{-2}$	$𝒦\left(A\right)$	1.15e2 (7)	8.76e2 (24)	6.94e3 (46)	5.54e4 (90)	4.43e5 (175)	3.54e6 (346)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (21)	2.44e1 (35)	9.57e1 (62)	3.81e2 (115)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.45 (11)	1.41 (10)	1.39 (9)	1.37 (8)	1.28 (6)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (8)	1.36 (10)	1.73 (13)	1.92 (14)	1.98 (14)	1.99 (13)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)	1.00 (4)
1	$𝒦\left(A\right)$	4.44 (6)	1.74e1 (21)	6.94e1 (40)	5.54e2 (78)	4.43e3 (153)	3.54e4 (302)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (21)	2.44e1 (35)	9.57e1 (62)	3.81e2 (115)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.45 (11)	1.41 (10)	1.39 (9)	1.33 (7)	1.28 (6)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (8)	1.36 (10)	1.73 (13)	1.92 (14)	1.98 (14)	1.99 (13)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)	1.00 (3)
${10}^2$	$𝒦\left(A\right)$	3.88e2 (6)	2.00e2 (22)	9.98e1 (44)	2.77e2 (80)	1.11e3 (143)	4.43e3 (273)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (21)	2.44e1 (35)	9.57e1 (62)	3.81e2 (115)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.45 (11)	1.41 (10)	1.39 (9)	1.33 (7)	1.28 (6)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (8)	1.36 (10)	1.73 (13)	1.92 (14)	1.98 (14)	1.99 (13)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)	1.00 (3)
${10}^4$	$𝒦\left(A\right)$	3.88e4 (6)	2.00e4 (22)	9.98e3 (47)	5.00e3 (89)	2.50e3 (163)	4.43e3 (295)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (21)	2.44e1 (35)	9.57e1 (62)	3.81e2 (115)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.45 (11)	1.41 (10)	1.39 (9)	1.33 (7)	1.28 (6)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (8)	1.36 (10)	1.73 (13)	1.92 (14)	1.98 (14)	1.99 (13)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)	1.00 (3)
${10}^6$	$𝒦\left(A\right)$	3.88e6 (9)	2.00e6 (22)	9.99e5 (51)	5.00e5 (96)	2.50e5 (180)	1.25e5 (331)
	$𝒦\left(B_{\mathrm{sgs}}A\right)$	1.10 (7)	2.11 (12)	6.54 (21)	2.44e1 (35)	9.57e1 (62)	3.81e2 (115)
	$𝒦\left(B_{\mathrm{fict}}A\right)$	1.32 (6)	1.45 (11)	1.41 (10)	1.39 (9)	1.33 (7)	1.28 (6)
	$𝒦\left(B_{\mathrm{add}}A\right)$	1.14 (8)	1.36 (10)	1.73 (13)	1.92 (14)	1.98 (14)	1.99 (13)
	$𝒦\left(B_{\mathrm{mul}}A\right)$	1.00 (5)	1.00 (5)	1.00 (5)	1.00 (4)	1.00 (4)	1.00 (3)

Equations

$$-\nabla \cdot (\kappa \nabla u)=f\text{mbox}in\Omega ,u{|}_{\partial \Omega }=0.$$

$$(u,v{)}_{0,\kappa ,G}:={\int }_G\kappa (\mathbf{x})uvd\mathbf{x},$$

$$a(u,v):={\int }_{\Omega }\kappa \nabla u\cdot \nabla vdx=(f,v),\forall v\in V.$$

$$V_h^K:=\{v\in H^1(K):v{|}_{\partial K}\in {\mathbb{B}}_k(\partial K),\Delta v\in {\mathbb{P}}_{k-2}(K)\},$$

$$V_h:=\{v\in V:v{|}_K\in V_h^K,\forall K\in {\mathcal{T}}_h\}.$$

$$\Vert \nabla v{\Vert }_{0,K}\le C{h}_K^{-1}\Vert v{\Vert }_{0,K},\forall v\in V_h^K.$$

$$(\nabla {\Pi }_h^{\nabla }v,\nabla p{)}_K=(\nabla v,\nabla p{)}_K,\forall p\in {\mathbb{P}}_k(K).$$

$$\sum _{i=1}^{n_v^K}v({\mathbf{x}}_i)$$

$${\int }_K{\Pi }_h^{\nabla }vd\mathbf{x}$$

$$a(u_h,v)=\sum _K{a}^K(u_h,v)=(f,v),\forall v\in V_h,$$

$$a^K(u_h,v)=a^K({\Pi }_h^{\nabla }{u}_h,{\Pi }_h^{\nabla }v)+a^K((I-{\Pi }_h^{\nabla })u_h,(I-{\Pi }_h^{\nabla })v).$$

$$a_h^K(u_h,v):=a^K({\Pi }_h^{\nabla }{u}_h,{\Pi }_h^{\nabla }v)+s^K((I-{\Pi }_h^{\nabla })u_h,(I-{\Pi }_h^{\nabla })v).$$

$$a_h^K(v,p)=a^K(v,p).$$

$$c_1{a}^K(v,v)\le a_h^K(v,v)\le c_2{a}^K(v,v),\forall v\in V_h^K.$$

$$a_h(u_h,v_h)=⟨f_h,v_h⟩,\forall v_h\in V_h,$$

$$(Av,w)=(v,w{)}_A:=a_h(v,w),\forall v,w\in V_h.$$

$$A{u}_h=b.$$

$$\mathcal{K}(A)\lesssim \mathcal{J}(\kappa )h^{-2},$$

$$(Av,v)$$

$$\le \sum _{K\in {\mathcal{T}}_h}{c}_2{a}^K(v,v)\le c_2\underset{\mathbf{x}\in \Omega }{\max }\kappa (\mathbf{x})\sum _{K\in {\mathcal{T}}_h}\Vert \nabla v{\Vert }_{0,K}^2$$

$$\lesssim h^{-2}\underset{\mathbf{x}\in \Omega }{\max }\kappa (\mathbf{x})(v,v).$$

$$\underset{\mathbf{x}\in \Omega }{\min }\kappa (\mathbf{x})(v,v)$$

$$\lesssim \sum _{K\in {\mathcal{T}}_h}{a}^K(v,v)$$

$$\lesssim (Av,v),\forall v\in V_h.$$

$$\mathcal{K}(A)={\lambda }_{\max }(A)/{\lambda }_{\min }(A)\lesssim h^{-2}\mathcal{J}(\kappa ).$$

$$a(u_h^c,v_h)=(f,v_h),\forall v_h\in V_h^c.$$

$$(A_{c}v,w)=(v,w{)}_{A_c}:=a(v,w),\forall v,w\in V_h^c.$$

$$({\Pi }^{t}v,w)=(v,\Pi w),\text{for }v\in V_h\text{ and }w\in V_h^c.$$

$$B_{\mathrm{add}}=\mathcal{R}+\Pi A_c^{-1}{\Pi }^t,$$

$$I-B_{\mathrm{mul}}A=(I-\mathcal{R}A)(I-\Pi A_c^{-1}{\Pi }^t)(I-\mathcal{R}A).$$

$$a_h(v,v)\le c_{0}s(v,v),\forall v\in V_h,$$

$$a_h(\Pi w,\Pi w)\le c_{1}a(w,w),\forall w\in V_h^c.$$

$$s(v_0,v_0)+a(w,w)\le c_2{a}_h(v,v).$$

Full text

YZ was supported in part by NSF DMS 1319110.

1. Introduction

In this paper, we present the auxiliary space preconditioning techniques ([20]) for the linear system arising from virtual element methods (VEMs) discretization on polytopal mesh for the following second order elliptic problems with jump coefficients:

[TABLE]

Here, $\Omega\subset{\mathbb{R}}^{d}$ ($d=2,3$) is an open polygonal domain and $f\in L^{2}(\Omega)$. We assume that the diffusion coefficient $\kappa\in L^{\infty}(\Omega)$ is piecewise constant with respect to an initial polytopal partition ${\mathcal{T}}_{0}$ of $\Omega$ but may have large jumps across the interface of the partition.

Recently introduced in [3], the VEMs have been widely used in the numerical approximation of various partial differential equations (PDEs) (see for example [1, 5, 9, 14]). The VEMs are generalizations of the classical finite element methods (FEMs), which permit the use of general polygonal and polyhedral meshes. Using polytopal meshes allows for more flexibility in dealing with complex computational domains or interfaces (cf. [13]). Some error estimates of the VEM are recently given in [8, 11]. For practical implementation of the the VEM for second order elliptic equations, we refer to [4, 17].

On the other hand, little work has been devoted to design efficient solvers for the VEM. In [6], the authors developed and analyzed balancing domain decomposition preconditioners for elliptic problems discretized by the virtual element method. In [2], the authors presented a $p$-version multigrid algorithm for the VEM discretization of 2D Poisson problem. The coarse spaces are constructed by decreasing progressively the polynomial degree of the VEM space. It was shown that the multigrid algorithm converges uniformly.

In order to design an efficient, robust and easily implementable preconditioner for VEM on polytopal mesh, we use the general auxiliary space preconditioning framework developed in [20]. This framework allows us to construct preconditioners in the lack of the hierarchical subspaces that is required in the classical (geometric) multilevel methods. The prominent examples of this framework include [15] for the Maxwell’s equations and [12] for the weak Galerkin method on simplicial mesh. The preconditioners we will develop here can be understood as two-level algorithms. The “fine” level problem is the VEM discretization on general polytopal mesh, and the “coarse” level problem is the standard conforming ${\mathbb{P}}_{1}$ finite element space on an auxiliary simplicial mesh. It is natural to choose the standard ${\mathbb{P}}_{1}$ finite element spaces as the coarse space for a couple of reasons: (1) the degrees of freedom in the coarse space are included in the VEM space – so asymptotically, the solutions on the coarse and fine spaces should be close to each other; (2) there is a lot of literature on developing efficient (and robust) solvers for the standard conforming ${\mathbb{P}}_{1}$ finite element discretization, so we can use any existing solvers/preconditioners as a coarse solver.

In this paper, we show that the auxiliary space preconditioners are robust with respect to both mesh size and the jumps in the coefficient. In addition, we show that the fictitious space preconditioner using the conforming ${\mathbb{P}}_{1}$ finite element space as the fictitious space also provides an efficient and robust preconditioner for the linear VEM discretization. We also remark that these preconditioners are very easy to implement. In particular, the matrix form of the intergrid transfer operator between the linear VEM and the coarse space is just the identity matrix.

The rest of this paper is organized as follows. In Section 2, we give basic notation and the virtual element discretization. Then in Section 3, we present the auxiliary space preconditioners and discuss its convergence. In Section 4, we present several numerical experiments in both 2D and 3D to verify the theoretical results. Finally, in Section 5 we summarize the results and discuss some possible extensions of these results. Throughout the paper, we will use the notation $x_{1}\lesssim y_{1}$, and $x_{2}\gtrsim y_{2}$, whenever there exist constants $C_{1},C_{2}$ independent of the mesh size $h$ and the coefficient $\kappa$ or other parameters that $x_{1}$, $x_{2}$, $y_{1}$ and $y_{2}$ may depend on, and such that $x_{1}\leq C_{1}y_{1}$ and $x_{2}\geq C_{2}y_{2}$.

2. Virtual Element Methods

Given any subset $G\subset{\mathbb{R}}^{d}$, we use standard notation $L^{2}(G)$ for the set of square integrable functions, with the inner product $(\cdot,\cdot)$ and the induced norm $\|\cdot\|_{0,G}$. We use the standard definition and notation of Sobolev spaces $H^{1}(G)$ and $H_{0}^{1}(G)$. The associated inner product, norm, and semi-norm in $H^{1}(G)$ are denoted by $(\cdot,\cdot)_{1,G}$, $\|\cdot\|_{1,G}$, and $|\cdot|_{1,G}$, respectively. We will also use $(\cdot,\cdot)_{0,\kappa,G}$ to denote the weighted $L^{2}$-inner product

[TABLE]

and $\|\cdot\|_{0,\kappa,G}$ as the weighted $L^{2}$-norm. For simplicity, when $G=\Omega$, we omit it from the norms (or the inner products). The variational formulation of the equation (1.1) reads: find $u\in V:=H_{0}^{1}(\Omega)$ such that

[TABLE]

Let ${\mathcal{T}}_{h}$ be a family of partitions of $\Omega$ into non-overlapping simple polytopal elements $K$. Here the subscript $h=h(K)$ denotes the diameter of the element $K$. Following [11], we make the following assumption on the polytopal mesh:

1. (A)

Each polytopal element $K\in{\mathcal{T}}_{h}$ has a “virtual triangulation” ${\mathcal{T}}_{K}$ of $K$ such that ${\mathcal{T}}_{K}$ is uniformly shape regular and quasi-uniform. Each edge of $K$ is an edge of certain element in ${\mathcal{T}}_{K}$.

Let $k\geq 1$ be an integer, and ${\mathbb{P}}_{l}(D)$ be the polynomial of degree less than or equal to $l$ on the domain $D$. For convenience, we denote ${\mathbb{P}}_{-1}(D)=\{0\}.$ On each polytopal element $K\in{\mathcal{T}}_{h}$, we define the local virtual finite element space:

[TABLE]

where ${\mathbb{B}}_{k}(\partial K):=\{v\in C^{0}(\partial K):v|_{e}\in{\mathbb{P}}_{k}(e),\;\forall e\subset\partial K\}.$ It is clear that $V_{h}^{K}\supset{\mathbb{P}}_{k}(K)$, and it may implicitly contain some other non-polynomial functions. The global virtual element space $V_{h}$ is then defined as:

[TABLE]

In the 2D case, any function $v\in V_{h}$ can be uniquely determined by the following degrees of freedom:

• •

the function values of $v$ at the vertices of ${\mathcal{T}}_{h}$;

• •

the function values of $v$ at the $k-1$ internal points of the $(k+1)$ Gauss-Lobatto quadrature rule on each edge $e$;

• •

the moments up to order $k-2$ of $v$ on each element $K\in{\mathcal{T}}_{h}$.

For a more detailed construction of the $V_{h}$ in both 2D and 3D, we refer to [4]. The following inverse inequality was shown in [11, Theorem 3.6] (see also [8, Lemma 2.19]).

Lemma 2.1**.**

There exists a constant $C$ depending only on the shape regularity and quasi-uniformity of ${\mathcal{T}}_{K}$ such that

[TABLE]

The VEM discretization is designed in such as way that the “virtual” basis functions are not necessary to be computed explicitly. Instead, they are understood through the degrees of freedom. Although the “virtual” functions $v\in V_{h}^{K}$ are defined implicitly, we can compute their Galerkin projection $\Pi_{h}^{\nabla}:V_{h}^{K}\to{\mathbb{P}}_{k}(K)$ onto the polynomial space ${\mathbb{P}}_{k}(K)$. That is,

[TABLE]

$\Pi_{h}^{\nabla}v$ can be uniquely determined by an additional constraint, e.g., by setting

[TABLE]

where $n_{v}^{K}$ is the total number of vertices of $K$, and $\boldsymbol{x}_{i}$ ($i=1,\cdots,n_{v}^{K}$) are the coordinates of the vertices.

Like the standard finite element methods, the VEM approximation of the equation (2.1) is to find $u_{h}\in V_{h}$ such that

[TABLE]

where $a^{K}(\cdot,\cdot):=a(\cdot,\cdot)|_{K}$ is the restriction of the bilinear form on $K$. For any $v\in V_{h}^{K}$, we may decompose $v=\Pi_{h}^{\nabla}v+(I-\Pi_{h}^{\nabla})v.$ By the orthogonality of $\Pi_{h}^{\nabla}$ and $(I-\Pi_{h}^{\nabla})$ in the energy norm, we have

[TABLE]

Now the first term can be evaluated exactly. For the second term, we replace it by some computable bilinear form $a^{K}((I-\Pi_{h}^{\nabla})u_{h},(I-\Pi_{h}^{\nabla})v)\approx s^{K}((I-\Pi_{h}^{\nabla})u_{h},(I-\Pi_{h}^{\nabla})v)$, to get

[TABLE]

The first term in the definition of $a_{h}^{K}(\cdot,\cdot)$ is called consistency term, and the second term is called the stabilizing term. For the well-posedness and convergence of the VEM, it is important that $a_{h}^{K}(\cdot,\cdot)$ satisfies

• •

Polynomial consistency: for any $v\in V_{h}^{K}$ and $p\in{\mathbb{P}}_{1}(K)$, it satisfies

[TABLE]

• •

Stability: There exist constants $c_{1}$ and $c_{2}$ independent of $\kappa$ and $h$ such that

[TABLE]

There are different choices for the bilinear form $s^{K}(\cdot,\cdot)$ that satisfy the polynomial consistency (2.2) and stability (2.3) conditions, see for example [3, 5]. In the numerical test performed in Section 4, we use the standard choice of $s^{K}$ as in [3] with weights $\kappa_{K}=\kappa(x)|_{K}$ on each element $K\in{\mathcal{T}}_{h}.$

Let $a_{h}(\cdot,\cdot):=\sum_{K\in{\mathcal{T}}_{h}}a_{h}^{K}(\cdot,\cdot).$ Then the VEM discretization of (2.1) reads: find $u_{h}\in V_{h}$ such that

[TABLE]

where $f_{h}$ is the $L^{2}$ projection of $f$ on the piecewise ${\mathbb{P}}_{k-2}(K)$ for each element $K\in{\mathcal{T}}_{h}.$ Let $A$ be the operator induced by the bilinear form $a_{h}(\cdot,\cdot)$, that is,

[TABLE]

Then solving (2.4) is equivalent to solve the linear system

[TABLE]

It is clear that the operator $A$ is symmetric and positive definite. Now we give an estimate of the condition number of $A$.

Lemma 2.2**.**

The condition number of the operator $A$, denoted by ${\mathcal{K}}(A)$, satisfies

[TABLE]

where ${\mathcal{J}}(\kappa):=\max_{\boldsymbol{x}\in\Omega}\kappa(x)/\min_{\boldsymbol{x}\in\Omega}\kappa(x)$ is the jump in the coefficient.

Proof.

By the stability of $a_{h}^{K}$ and inverse inequality Lemma 2.1, we have

[TABLE]

Since $A$ is symmetric and positive definite with respect to $(\cdot,\cdot)$, the above inequality implies that $\lambda_{\max}(A)\lesssim h^{-2}\max_{\boldsymbol{x}\in\Omega}\kappa(\boldsymbol{x})$.

On the other hand, by the Poincaré inequality and the stability of $a_{h}^{K}$, we obtain

[TABLE]

This implies $\lambda_{\min}(A)\gtrsim 1/\min_{\boldsymbol{x}\in\Omega}\kappa(\boldsymbol{x})$. Therefore

[TABLE]

This completes the proof. ∎

This lemma implies that the linear algebraic system resulting from the VEM discretization (2.4) is ill-conditioned with the condition number depends on both the mesh size and the jump in the coefficient. This dependence can be observed from the numerical experiments presented in Section 4. Thus the linear system (2.5) is difficult to solve using the classic iterative methods such as Jacobi, Gauss-Seidel or CG, without effective preconditioners. The main purpose in this paper is the develop efficient auxiliary space preconditioners for the linear system (2.5).

Remark 2.3**.**

Based on the norm equivalences obtained in [11], the results in this paper can be applied to the other type of VEM discretizations directly (namely, with different form of stabilization term).

3. Auxiliary Space Preconditioner

In this section, we construct an auxiliary space preconditioner for solving the discrete system of equation (2.5). To do that, we need to introduce an auxiliary space. For each $K\in{\mathcal{T}}_{h}$, we introduce an auxiliary triangulation such that each edge of $K$ is an edge of some element in this triangulation. This can be done using the Delaunay triangulation. This leads to a conforming triangulation of the whole domain $\Omega.$ Moreover, with the Assumption (A) on ${\mathcal{T}}_{h}$, the resulting triangulation is quasi-uniform. On this triangulation, we define a standard conforming ${\mathbb{P}}_{1}$ finite element subspace $V^{c}_{h}$. We consider the auxiliary problem: find $u^{c}_{h}\in V^{c}_{h}$ such that

[TABLE]

Similiary, let $A_{c}$ be the operator induced by the bilinear form $a(\cdot,\cdot)$, that is,

[TABLE]

The operator $A_{c}$ is symmetric and positive definite.

To construct the auxiliary space preconditioner, we use $V_{h}$ as the “fine” space and $V_{h}^{c}$ as the “coarse” space. Since $A_{c}$ is the conforming piecewise linear finite element discretization of the equation (2.1), the “coarse” problem in $V_{h}^{c}$ can be solved by many existing efficient solvers such as the standard multigrid, domain decomposition or algebraic multigrid (AMG) solvers (see for example [21, 22]). Next, on the fine space, we define a “smoother” ${\mathcal{R}}:\>V_{h}\rightarrow V_{h}$, which is symmetric positive definite. For example, ${\mathcal{R}}$ could be a Jacobi or symmetric Gauss-Seidel smoother. We denote $s(\cdot,\cdot)$ as the bilinear form corresponding to ${\mathcal{R}}^{-1}$. Finally, to connect the “coarse” space $V_{h}^{c}$ with the “fine” space $V_{h}$, we need a “prolongation” operator $\Pi:V_{h}^{c}\to V_{h}$. A “restriction” operator $\Pi^{t}:V_{h}\rightarrow V_{h}^{c}$ is then defined as

[TABLE]

Then, the auxiliary space preconditioner $B:\>V_{h}\rightarrow V_{h}$ is given by

[TABLE]

With the notation introduced above, we have the following theorem for the auxiliary space preconditioners.

Theorem 3.1** (cf. [20, 15]).**

Let the auxiliary space preconditioner defined above satisfy the following conditions:

• (C0)

The smoother ${\mathcal{R}}$ satisfies that

[TABLE]

where $c_{0}>0$ is independent $h$ and $\kappa$.

• (C1)

The operator $\Pi:V_{h}^{c}\to V_{h}$ is stable in the sense that there exists a constant $c_{1}>0$ independent of $h$ and $\kappa$ such that

[TABLE]

• (C2)

For any $v\in V_{h}$, there exist $w\in V_{h}^{c}$ and $v_{0}\in V_{h}$ such that $v=v_{0}+\Pi w$ such that

[TABLE]

Then

[TABLE]

Remark 3.2**.**

1. (1)

In the definition of the auxiliary space preconditioners (3.2)-(3.3), the coarse space solver $A_{c}^{-1}$ could be replaced by an inexact solver ${\mathcal{B}}_{c}\approx A_{c}^{-1}$ (see for example **[15, Corollary 2.3]**). For simplicity, we use the exact solver in the analysis. 2. (2)

In  (3.2), if we ignore the smoother ${\mathcal{R}}$, then the preconditioner is usually called fictitious space preconditioner (cf. **[16]**). In the sequel, we denote

[TABLE]

The auxiliary space preconditioner can be viewed as a generalization of the fictitious space preconditioner by a special choice of the “fictitious space”. In particular, the fictitious space is defined as a product space having $V_{h}$ itself as one of the component. This makes it easier to construct the map from the fictitious space to the original space, which is required to be surjective.

The proof of Theorem 3.1 amounts to verifying three conditions (C0)-(C2). The smoothing property (C0) is given by the following lemma. The proof of the smoothing properties are standard. The proof can be carried out by showing the estimates for the Jacobi smoother, as the Jacobi smoother and the symmetric Gauss-Seidel smoother are spectrally equivalent for any symmetric positive definition matrix (see for example [19, Proposition 6.12] or [23, Lemma 3.3].) We omit the details here.

Lemma 3.3**.**

Let $s(\cdot,\cdot)$ be the bilinear form associated to Jacobi, or symmetric Gauss-Seidel smoother. Then for any $v\in V_{h}$, we have the following estimates

[TABLE]

Here, the constants are independent of coefficient and mesh size.

For the intergrid transfer operator $\Pi:V_{h}^{c}\to V_{h},$ we defined on each element $K\in{\mathcal{T}}_{h}$

[TABLE]

for each $w\in V_{h}^{c}$; namely, $\Pi w$ is the harmonic extension of $w|_{\partial K}\in{\mathbb{B}}_{1}(\partial K).$ We remark that when $k=1$, it is not necessary to solve (3.6) in order to construct the operator $\Pi.$ In fact, in this case the matrix representation of $\Pi$ is just the identity matrix due to the corresponding degrees of freedom for $V_{h}$ and $V_{h}^{c}$. This intergrid transfer operator defined in (3.6) satisfies the following properties.

Lemma 3.4**.**

For any $w\in V_{h}^{c}$, the operator $\Pi:V_{h}^{c}\to V_{h}$ defined by (3.6) satisfies the following properties:

[TABLE]

Proof.

By the stability assumption (2.3) we only need to show that

[TABLE]

By the definition of $\Pi w$ and trace inequality, we have

[TABLE]

This implies (3.7).

To prove (3.8), we can view $w\in V_{h}^{c}$ as piecewise linear interpolation of $\Pi w$ on each element $K\in{\mathcal{T}}_{h}$. Therefore, by the stability of $\Pi$ (3.7) we have

[TABLE]

The estimate (3.8) then follows. ∎

Inequality (3.7) is exactly the condition (C1) in Theorem 3.1. It remains to verify the stable decomposition condition (C2). This is achieved by the following lemma.

Lemma 3.5**.**

There exists a linear operator $P:V_{h}\to V_{h}^{c}$ such that

[TABLE]

Proof.

Let $P=I_{K}:V_{h}^{K}\to V_{h}^{c}(K)$ be the local interpolation operator on the subdomain $K\in{\mathcal{T}}_{h}$ such that $I_{K}(v)(\boldsymbol{x}_{i})=v(\boldsymbol{x}_{i})$ for each $v\in V_{h}^{K}$ and each vertex $\boldsymbol{x}_{i}\in K$. Note this local interpolation operator is not well-defined for general $H^{1}(K)$, but it is well-defined on the VEM space $V_{h}$. Then it is obvious that

[TABLE]

and hence $P$ satisfies the stability (3.9).

To prove (3.10), by triangle inequality

[TABLE]

where in the second inequality we used the approximation (3.11) of $P$ for the first term, and we used (3.8) for the second term. The estimate (3.10) then follows by (3.9). This completes the proof. ∎

Lemma 3.5 implies that for any $v\in V_{h}$, we can decompose it as $v=v_{0}+\Pi w,$ where $w=Pv\in V_{h}^{c}$ and $v_{0}=v-\Pi w\in V_{h}$. With the smoothing property (3.5), we get

[TABLE]

This verifies the condition (C2). In summary, by the auxiliary space preconditioner framework Theorem 3.1, we have the following main theorem.

Theorem 3.6**.**

The auxiliary space preconditioner (both additive and multiplicative) defined in (3.2) (3.3) satisfies:

[TABLE]

where the constant $C$ is independent of the mesh size $h$ and $\kappa$.

Note that in general, the intergrid operator defined by (3.6) is not surjective. In fact, there is no surjective intergrid operator from $V_{h}^{c}$ to $V_{h}$ when $k\geq 2$. Hence, when $k\geq 2$ it is necessary to use the auxiliary space preconditioner. However, when $k=1$ we have the following corollary for the fictitious space preconditioner $B_{{\rm fict}}$.

Corollary 3.7**.**

When $k=1$, namely the linear VEM, the fictitious space preconditioner $B_{{\rm fict}}=\Pi A_{c}^{-1}\Pi^{t}$ satisfies:

[TABLE]

where the constant $C$ is independent of the mesh size and $\kappa$.

Proof.

When $k=1$, the operator $\Pi$ defined in (3.6) is surjective and satisfies (3.7) and (3.9). Therefore, the conclusion follows immediately by the fictitious space lemma (cf. [15, Theorem 2.2]). ∎

4. Numerical Experiments

In this section, we present several numerical experiments in both 2D and 3D to verify the result in Theorem 3.6 on the performance of the proposed preconditioners. In all these tests, we use 2-sweeps symmetric Gauss-Seidel smoother. The stopping criteria is $\|r_{k}\|/\|r_{0}\|<10^{-12}$ for the PCG algorithm, where $r_{k}=f-Au_{k}$ is the residual. For the coarse solver, we use the AMG algorithm implemented in $i$FEM [10].

4.1. 2D Examples

In the first example, we consider the model problem (1.1) in the unit square $\Omega=[0,1]^{2}$ with constant coefficient $\kappa=1$. Figure 2 is an example of the polytopal mesh of the unit square domain (with 100 elements) generated using yMesher [18], and Figure 2 is the corresponding Delaunay triangular mesh. The VEM discretization is defined on the polytopal mesh (cf. Figure 2), while the auxiliary space using the standard conforming ${\mathbb{P}}_{1}$ finite element discretization is defined on the corresponding triangular mesh (cf. Figure 2).

Tables 1 shows the estimated condition numbers (the number of PCG iterations) for the additive and multiplicative preconditioned systems.

For comparison, we also include the estimated condition numbers ${\mathcal{K}}(A)$, ${\mathcal{K}}(B_{{\rm sgs}}A)$ and ${\mathcal{K}}(B_{{\rm fict}}A)$, where $B_{{\rm sgs}}$ is the (2-sweep) symmetric Gauss-Seidel preconditioner (same below) and $B_{{\rm fict}}$ is the fictitious space preconditioner using the conforming FEM. As we can observe from this table, while the condition numbers ${\mathcal{K}}(A)$, ${\mathcal{K}}(B_{{\rm sgs}}A)$ increase as the mesh is refined. The condition number ${\mathcal{K}}(B_{{\rm fict}}A)$ increase slightly as the mesh is refined. On the other hand, the condition numbers of $\mathcal{K}(B_{{\rm add}}A)$ and $\mathcal{K}(B_{{\rm mul}}A)$ are uniformly bounded.

In the second test, we consider the problem with jump coefficients. The coefficients $\kappa$ are generated randomly on each polygon element (see Figure 3 for an example of the coefficient distribution with 100 elements, the integer in each polygonal element is the magnitude of the coefficient.).

Note that the coefficient settings are different in different polytopal meshes. Tables 2 shows the estimated condition numbers (the number of PCG iterations). Here, “-” means the PCG algorithm failed to converge after 1200 iterations. As we can see from this table, while ${\mathcal{K}}(A)$ and ${\mathcal{K}}(B_{{\rm sgs}}A)$ increase dramatically, the condition number $\mathcal{K}(B_{{\rm add}}A)$ and $\mathcal{K}(B_{{\rm mul}}A)$ are nearly uniformly bounded. These observations verify the estimate given in Theorem 3.6.

In the third test, we consider the performance of the preconditioners for voronoi meshes which violate the assumption (A) (see Figure 4 for an example of 100 polygons). As we can observe from this figure, the aspect ratio for some polygons are quite high – thus the partition is no longer quasi-uniform. Similar to before, we use Delaunay triangulation of this mesh to construct the auxiliary space.

Table 3 shows the estimated condition numbers with the number of PCG iterations for different preconditioners. As we can see from this table, both additive and multiplicative auxiliary space preconditioners are still robust with respect to the problem size.

4.2. 3D Example

Now we consider the model problem (1.1) in a 3D cubic domain $\Omega=[0,1]^{3}$. We subdivide the domain into hexagonal elements (cubes) with mesh size $h$ at each level. The VEM discretization is defined on the hexagon mesh. For the auxiliary space, we further divide each hexagon into six tetrahedrons to construct the auxiliary mesh and to define the ${\mathbb{P}}_{1}$ conforming finite element discretization on this auxiliary mesh (see for example, Figure 5).

In this example, we test various discontinuous coefficient settings. Let $\Omega_{1}=[0.25,0.5]^{3}$ and $\Omega_{2}=[0.5,0.75]^{3}$ (see Figure 5). We set the coefficient $\kappa|_{\Omega_{1}\cup\Omega_{2}}=\kappa_{1}=10^{k}$ (with $k=-6,-4,-2,0,2,4,6$) and $\kappa|_{\Omega\setminus(\Omega_{1}\cup\Omega_{2})}=1$.

Table 4 presents the estimated condition number of the preconditioned systems with respect to different choice of $\kappa_{1}$ and mesh size.

As we can see from Table 4, the condition number of $A$ depends on both the coefficient $\kappa$ and the mesh size. On the other hand, both the fictitious space preconditioner and the auxiliary space preconditioners (additive or multiplicative) are efficient and robust with respect to jumps in the coefficient $\kappa$ and the mesh size. These results justify Theorem 3.6 and Corollary 3.7.

5. Conclusion

In this paper, we developed the additive and multiplicative auxiliary space preconditioners for solving the linear system arising from the virtual element methods discretization on polytopal meshes for the second order elliptic equation with jump coefficients. We used an auxiliary simplicial triangulation to construct the coarse space. The auxiliary space preconditioners consist of a smoother and a coarse space correction. We showed that the condition numbers of the preconditioned systems are uniformly bounded, independent of the problem size and the jump in coefficients. For the linear VEM discretization, we also showed that the fictitious space preconditioner is also a robust and efficient preconditioner. Numerical experiments were presented to demonstrate the performance of these preconditioners.

In the analysis, we assumed that the polytopal mesh should satisfy (A), that is, there exist a quasi-uniform auxiliary triangulation. However, the numerical experiments demonstrated that these preconditioners also perform uniformly for voronoi meshes which violate the quasi-uniform assumption. This indicates that we might be able to relax the quasi-uniformity condition. This is a topic of ongoing research.

Acknowledgement

This work is supported in part by NSF DMS 1319110.