The Discontinuous Galerkin (DG) method was introduced in 1973 by Reed and Hill for the discretization of hyperbolic equations ReedHill73 . Extensions of the method were quickly proposed to deal with elliptic and parabolic problems: some of the most relevant works include Arnold Ar1982 , Baker Baker77 , Nitsche Nitsche and Wheeler Wh1978 , whose contributions put the basis for the development of the interior penalty DG methods. In the last 40 years the scientific and industrial community has shown an exponentially growing interest in DG methods - see for example CockKarnShu00 ; DiPiEr ; HestWar ; Riviere for an overview. On one side, the features of DG methods have been naturally enhanced by the recent development of High Performance Computing technologies as well as the growing request for high-order accuracy. In particular, as the discrete polynomial space can be defined locally on each element of the mesh, DG methods feature a high-level of intrinsic parallelism. Moreover, the local conservation properties and the possibility to use meshes with hanging nodes make DG methods interesting also from a practical point of view. Recently, it has been shown that DG methods can be extended to computational grids characterized by polytopic elements, cf. Ref. AnBrMa2008 ; AnFaRuVe2016 ; AntGiaHou13 ; AntGiaHou14 ; AnHoHuSaVe2017 ; AnMa17 ; BaBoCoSu14 ; Basetal12 ; BaBoetal12 ; Canetal14 ; GiHo2014 ; LiDaVaYo2014 ; WiKuMiTaWeDa2013 . In particular, the efficient approach presented in Canetal14 is based on defining a local polynomial discrete space by making use of the bounding box of each element GiHo2014_2 : this technique together with a careful choice of the discontinuity penalization parameter permits the use of polytopal elements which can be characterized by faces of arbitrarily small measure and as shown in CaDoGe2016 , see also AnHoHuSaVe2017 , possibly by an unbounded number of faces.
On the other hand, the development of fast solvers and preconditioners for the linear system of equations arising from high-order DG discretization is been developed. A recent strand of the literature has focused on multilevel techniques, including Schwarz domain decomposition methods, cf. Ref. AntGiaHou14 ; AnHoSm16 , and two-level and multigrid techniques, cf. Ref. AnHoHuSaVe2017 ; AntSarVer15 . The efficiency of those methods is more evident in the case of polygonal grids, because the flexibility of the element shape couples very well with the possibility to easily define agglomerated meshes, which is the key ingredient for the developing of multigrid algorithms. In AnHoHuSaVe2017 a two-level scheme and W-cycle multigrid method is developed to solve the linear system of equations arising from high-order discretization introduced in Canetal14 . One iteration of the proposed methods consists of an iterative application of the smoothing Richardson operator and the subspace correction step. In particular, the latter is based on a nested sequence of discrete polynomial spaces where the underlying polytopal grid of each subspace is defined by agglomeration. While being faster than other classical iterative methods, the agglomeration approach presents itself some limitations. When the finest grid is unstructured and characterized by polytopic elements, there is the possibility that its very small edges could be inherited by the coarser levels until the one where the linear system is solved with a direct method. In this case the presence of small faces negatively affects the condition number of the associated matrix: indeed, according to Canetal14 , the discontinuity penalization parameter is defined locally in each face as the inverse of its measure.
In this paper we aim to overcome this issue by solving the same linear system through a multilevel method characterized by a sequence of non-nested agglomerated meshes in order to make sure that the number of faces of the agglomerates does not blows up as the number of levels of our multigrid method increases. This can be achieved for example based on employing edge-coarsening techniques in the agglomeration procedures. The flexibility in the choice of the computational sub-grids leads to the definition of a non-nested multigrid method characterized by a sequence of non-nested multilevel discrete spaces, cf. Ref. BrVe1990 ; Zh1990 ; ZhZh1997 , and where the discrete bilinear forms are chosen differently on each level, cf. Ref. GoKa2003 ; GoKa2003_2 ; MoZh1995 . The first non-nested multilevel method was introduced by Bank and Dupont in BaDu1981 ; a generalized framework was developed by Bramble, Pasciak and Xu in BrPaXu1991 , and then widely used in the analysis of non-nested multigrid iterations, cf. Ref. Br1993 ; BrKwPa1994 ; BrPa1992 ; BrPa1993 ; BrPa1994 ; BrZh2001 ; GoPa2000 ; ScZh1992 ; XuCh2001 ; XuLiCh2002 . The method of BrPaXu1991 , to whom we will refer as the BPX multigrid framework, is able to generalize also the multigrid framework that we will develop in this paper, but the convergence analysis relies on the assumption that , which might not be guaranteed in the DG setting, as we will see in Sect. 4.2. Here and are two bilinear forms suitably defined on two consecutive levels, and is the prolongation operator whose definition is not trivial, differently from the nested case. For this reason the convergence analysis will be presented based on employing the abstract setting proposed by Duan, Gao, Tan and Zhang in DuanGaoTanZhang , which permits to develop a full analysis of V-cycle multigrid methods in a non-nested framework relaxing the hypothesis . We will prove that our V-cycle scheme with non-nested spaces converges uniformly with respect to the discretization parameters provided that the number of smoothing steps, which depends on the polynomial approximation degree , is chosen sufficiently large. This result extends the theory of AnHoHuSaVe2017 where W-cycle multigrid methods for high-order DG methods with nested spaces where proposed and analyzed.
The paper is organized as follows. In Sect. 2 we introduce the interior penalty DG scheme for the discretization of second-order elliptic problems on general meshes consisting of polygonal/polyhedral elements. In Sect. 3, we recall some preliminary analytical results concerning this class of schemes. In Sect. 4 we define the multilevel BPX framework for the V-cycle multigrid solver based on non-nested grids, and present the convergence analysis of the algorithm. The main theoretical results are validated through a series of numerical experiments in Sect. 5. In Sect. 6 we propose an improved version of the algorithm, obtained by choosing a smoothing operator based on a domain decomposition preconditioner.
2 Model problem and its DG discretization
We consider the weak formulation of the Poisson problem, subject to a homogeneous Dirichlet boundary condition: find such that
with , , a convex polygonal/polyhedral domain with Lipschitz boundary and . The unique solution of problem (1) satisfies
In view of the forthcoming multigrid analysis, let be a sequence of tessellation of the domain , each of which is characterized by disjoint open polytopal elements of diameter , such that , . The mesh size of is denoted by . To each we associate the corresponding discontinuous finite element space , defined as
where denotes the local space of polynomials of total degree at most on .
For the sake of brevity we use the notation to mean , where is a constant independent from the discretization parameters. Similarly we write in lieu of , while is used if both and hold.
A suitable choice of and leads to the -multigrid non-nested schemes. This method is based on employing, from one side, a set of non-nested partitions , such that the coarse level is independent from , with the only constrain
from the other side we assume that the polynomial degree vary from one level to another such that
Additional assumptions on the grids are outlined in the following paragraph.
2.1 Grid assumptions
For any , we define the faces of the mesh [T][j], , as the intersection of the -dimensional facets of neighbouring elements. This implies that, for , a face always consists of a line segment, however for , the faces of [T][j] are general shaped polygons. Thereby, we assume that each facets of an element may be subdivided into a set of co-planar -dimensional simplices and we refer to them as faces. In order to introduce the DG formulation, it is helpful to distinguish between boundary and interior element faces, denoted as and , respectively. In particular, we observe that for , while for any we assume that , where are two adjacent elements in . Furthermore, we denoted as the set of all mesh faces of . With this notation, we assume that the sub-tessellation of element interfaces into -dimensional simplices is given. Moreover, assume that the following assumptions hold, cf. CaDoGe2016 ; Canetal16 .
For any , given there exists a set of non-overlapping d-dimensional simplices , , such that for any face it holds that for some l, it holds , and the diameter of can be bounded by
For any , we assume that , where is the dimension of .
Every polytopic element , admits a sub-triangulation into at most shape-regular simplices , for some , such that and
Let , denote a covering of consisting of shape-regular d dimensional simplices . We assume that, for any , there exists such that and
are required for the inverse estimates of Lemma5 and Theorem 3.2. Assumption 2.4 guarantees the validity of the approximation result and error estimetes of Lemma 4 and Theorem 3.1, respectively.
2.2 DG formulation
In order to introduce the DG discretization of (1), we firstly need to define suitable jump and average operators across the faces , . Let and be sufficiently smooth functions. For each internal face , such that , let
be the outward unit normal vector to, and let and be the traces of the functions and on from , respectively. The jump and average operators across are then defined as follows:
cf. Arnetal01 . With this notation, the bilinear form corresponding to the symmetric interior penalty DG method on the -th level is defined by
where denotes the interior penalty stabilization function, which is defined by
with independent of , and , and is the lifting operator on the space , defined as
We refer to Arnetal01 for more details.
Here, the formulation with the lifting operators allows to introduce the discrete gradient operator , defined as
where is the piecewise gradient operator on the space . The role of will be clarified in Sect. 4.2.
The goal of this paper is to develop non-nested V-cycle multigrid schemes to solve the following problem posed on the finest level : find such that
By fixing a basis for , i.e. , formulation (16) results in the following linear system of equations
where is the vector of unknowns.
3 Preliminary results
In this section we recall some preliminary results which form the basis of the convergence analysis presented in the next section.
Assume that the sequence of meshes , satisfies Assumption 2.1 and let , then the following bound holds
where is the diameter of and is a positive number.
Assume that the sequence of meshes satisfies Assumption 2.1 and let . Then, the following bound holds
We refer to CaDoGe2016 for the proof.
On each discrete space , , we consider the following DG norm:
The well-posed of the DG formulation is established in the following lemma.
The following continuity and coercivity bounds, respectively, hold
Next, we recall the following approximation result, which is an analogous bound presented in (Canetal14, , Theorem 5.2).
Let Assumption 2.4 be satisfied, and let such that, for some , for each . Then there exists a projection operator such that
The result presented in Lemma 4 leads to the following error bounds for the underlying interior penalty DG scheme. The error in the energy norm has been proved in Canetal14 , see also CaDoGe2016 . -estimates can be found in AnHoHuSaVe2017 .
We point out that the bounds in Theorem 3.1 are optimal in and suboptimal in of a factor and for the -norm and the -norm, respectively. Optimal error estimates with respect to can be shown, for example, by using the projector of GeoSu03 for quadrilateral meshes providing the solution belongs to a suitable augmented Sobolev space. The issue of proving optimal estimates as the ones in GeoSu03 on polytopic meshes is an open problem and it is under investigation. In the following, we will write:
where , , and for optimal and suboptimal estimates, respectively.
We also need to introduce an appropriate inverse inequality, cf. Canetal16 .
Thanks to the inverse estimate of Lemma 5
, it is possible to obtain the following upper bound on the maximum eigenvalue of . We refer toAntHou for a similar result on standard grids, and to AnHoHuSaVe2017 for its extension to polygonal grids.
4 The BPX-framework for the V-cycle algorithms
The analysis presented in this section is based on the general multigrid theoretical framework already employed and developed in BrPaXu1991 for non-nested spaces and non-inherited bilinear forms. In order to develop a geometric multigrid, the discretization at each level follows the one already presented in AntSarVer15 , where a W-cycle multigrid method based on nested subspaces is considered. The key ingredient in the construction of our proposed multigrid schemes is the inter-grid transfer operators.
Firstly, we introduce the operators , defined as
and we denote as the maximum eigenvalue of . Moreover, let be the identity operator on level . The smoothing scheme, which is chosen to be the Richardson iteration, is then characterized by the following operators:
The prolongation operator connecting the coarser space to the finer space is denoted by . Since the two spaces are non-nested, i.e. , it cannot be chosen as the ”natural injection operator”. The most natural way to define the prolongation operator is the -projection, i.e.
The restriction operator is defined as the adjoint of with respect to the -inner product, i.e.,
For our analysis, we also need to introduce the operator such that:
where is defined as . Given an initial guess , and choosing parameters , the multigrid V-cycle iteration algorithm for the approximation of is outlined in Algorithm 1. In particular, represents the approximate solution obtained after one iteration of our non-nested V-cycle scheme, which is defined by induction: if we consider the general problem of finding such that
with and , then represents the approximate solution of (35) obtained after one iteration of the non-nested V-cycle scheme with initial guess and number of pre-smoothing and post-smoothing steps, respectively. The recursive procedure is outlined in Algorithm 2, where we also observe that on the level the problem is solved by using a direct method.
4.1 Convergence analysis
We first define the following norms on each discrete space
To analyze the convergence of the algorithm, for any we set and let be its adjoint respect to . Following DuanGaoTanZhang , we make three standard assumptions in order to prove the convergence of Algorithm 1:
Stability estimate: such that
Regularity-approximation property: such that
Smoothing property: such that
The convergence analysis of the V-cycle method is described by the following theorem that gives an estimate for the error propagation operator, which is defined as
We refer to DuanGaoTanZhang for the proof of Theorem 4.1 in an abstract setting. In the following, we prove the validity of Assumptions 1, 2 and 3 for the algorithm presented in this section. We start with a two-level approach, i.e. , so we will consider the two-level method for the solution of (16), based on two spaces . The generalization to the V-cycle method will be given at the end of this section.
4.2 Verification of Assumption 1
In order to verify Assumption 1 for the two-level method we first show a stability result of the prolongation operator . In the following, we also consider the -projection operator on the space defined as
From the definition of given in (31), it holds .
Moreover, we need the following approximation result which shows that any can be approximated by an -function, see AnHoSm16 . Let be the discrete gradient operator (15) introduced in Remark 4, and consider the following problem: , find such that
It is shown in AnHoSm16 that possesses good approximation properties in terms of providing an -conforming approximant of the discontinuous function :
Let be a bounded convex polygonal/polyhedral domain in , . Given , we write to be the approximation defined in (43). Then, the following approximation and stability results hold:
We make use of the previous result in order to show the following stability result of the prolongation operator:
There exists a positive constant , independent of the mesh size such that
Let , by the definition of the DG-norm (20), we need to estimate:
We next bound each of the two terms on the right hand side. For the first one let be defined as in (43). Then:
where we have added and subtracted the terms and
. The second term of the right hand above side can be estimated using the interpolation bounds of Lemma4, the Poincaré inequality for and the second bound of (44):
By adding and subtracting to we obtain
Using Lemma 4 and the Poincaré inequality we have
whereas the term can be estimate as follow: