Fast multigrid solvers for conforming and non-conforming multi-patch Isogeometric Analysis

# Fast multigrid solvers for conforming and non-conforming multi-patch Isogeometric Analysis

Stefan Takacs Johann Radon Institute for Computational and Applied Mathematics (RICAM),
###### Abstract

Isogeometric Analysis is a high-order discretization method for boundary value problems that uses a number of degrees of freedom which is as small as for a low-order method. Standard isogeometric discretizations require a global parameterization of the computational domain. In non-trivial cases, the domain is decomposed into patches having separate parameterizations and separate discretization spaces. If the discretization spaces agree on the interfaces between the patches, the coupling can be done in a conforming way. Otherwise, non-conforming discretizations (utilizing discontinuous Galerkin approaches) are required. The author and his coworkers have previously introduced multigrid solvers for Isogeometric Analysis for the conforming case. In the present paper, these results are extended to the non-conforming case. Moreover, it is shown that the multigrid solves get even more powerful if the proposed smoother is combined with a (standard) Gauss-Seidel smoother.

###### keywords:
Isogeometric Analysis, Multi-patch domains, Symmetric interior penalty discontinuous Galerkin
journal: arXiv\newdefinition

remarkRemark \newproofproofProof

## 1 Introduction

Isogeometric Analysis (IgA), see Hughes:2005 (), is a spline based approach for approximating the solution of a boundary value problem (BVP). One of the big strengths of IgA is that it has the approximation power of a high-order method while the number of degrees of freedom behaves basically like for a low-order method. To obtain this behavior, we have to be able to increase the spline degree while we simultaneously increase the smoothness. In IgA, this is typically called -refinement and leads to spline based discretizations.

In IgA, the computational domain (usually called physical domain) is parameterized using spline or NURBS functions. Since it might be too restrictive to parameterize the whole computational domain using just one global (smooth) geometry function, one typically represents the physical domain as the union of subdomains, in IgA called patches. Then, each of the patches is parameterized by its own geometry function (multi-patch IgA).

On each patch, a space of trial and test functions is introduced. The simplest approach to set up such a function space is to use tensor-product B-splines on the unit square or unit cube (parameter domain) and to use the geometry function to map them onto the physical domain or, in the multi-patch case, onto one of the patches. If we set up the function spaces such that the basis functions on the interfaces between the patches agree, we can use conforming discretizations. Approximation errors (cf. Hughes:2005 (); Bazilevs:BeiraoDaVeiga:Cottrell:Hughes:Sangalli:2006 (); BeiraoDaVeiga:Buffa:Rivas:Sangalli:2011 (); BeiraoDaVeiga:Buffa:Sangalli:Vazquez:2014 (); Takacs:Takacs:2015 (); Floater:Sande:2017 (); Sande:Manni:Speleers:2018 (); deLaRiva:Rodrigo:Gaspar:2018 (); Tielen:Moller:Goddeke:Vuik:2019 () and many others) and multigrid solvers (cf. Gahalaut:Kraus:Tomar:2013 (); Speleers:2015 (); Hofreither:Takacs:Zulehner:2017 (); Hofreither:Takacs:2017 () and others) for such discretizations have been previously discussed. Since we are interested in -refinement, we need results that are explicit in the spline degree. For the single-patch case, such error estimates have originally been given in Takacs:Takacs:2015 () and later improved in Sande:Manni:Speleers:2018 (). In Hofreither:Takacs:2017 (), a robust single-patch multigrid solver has been proposed and analyzed based on the error estimates from Takacs:Takacs:2015 (). In Takacs:2018 (), both the approximation error estimates and the multigrid solver have been extended to the conforming multi-patch case. These results are the foundation of the present paper.

If conforming discretizations are not feasible, discontinuous Galerkin (dG) approaches are possible. One standard dG approach is the Symmetric Interior Penalty discontinuous Galerkin (SIPG) method, see Arnold:1982 (); Arnold:Brezzi:Cockburn:Marini:2002 (). Already in LMMT:2015 (); LT:2015 (), it has been proposed to utilize these approaches to couple patches in IgA. Recently, also the dependence of the approximation error on the spline degree has been analyzed, see Takacs:2019 (). It was not possible to show that the approximation error is robust in the spline degree but it could be proven that it only grows logarithmically.

(Robust) multigrid solvers for such non-conforming discretizations are not known so far. In the present paper, it is shown how the multigrid solver from Takacs:2018 () can be extended to SIPG discretizations; we observe – as in Takacs:2018 () – that the numerical experiments show both robustness in the grid size and the spline degree. For completeness, we also show how to extend the convergence analysis from Takacs:2018 () to SIPG discretizations. It is worth noting that there are alternative solvers for multi-patch IgA, like FETI-type approaches, cf. KleissEtAl:2012 (); Hofer:2017 () and others, overlapping Schwarz type methods, cf. BeiraoDaVeiga:Cho:Pavarino:Scacchi:2012 (), or BDDC methods, cf. BeiraoDaVeiga:Cho:Pavarino:Scacchi:2013 (); most of them, however, have not been worked out for the non-conforming case.

Note that the idea behind the proposed subspace corrected mass smoother is that the boundary value problem on the physical domain (on one patch) can be well approximated by a boundary value problem on the parameter domain. Thus, the tensor-product structure on the parameter domain can be used. This is true if the geometry function is not too distorted. Otherwise, the convergence behavior suffers significantly. The same behavior can be observed by other fast solvers that are based on the same idea, cf. the fast diagonalization method Sangalli:2016 (). Here, the authors have improved their method by incorporating the geometry information into the preconditioner, cf. Montardini:Sangalli:Tani:2018 (). For the multigrid setting, it has turned out that one can overcome these problems quite well if the subspace corrected mass smother is combined with a Gauss-Seidel smoother (hybrid smoother) since both approaches have strengths that seem to be somewhat orthogonal to each other (robustness in spline degree vs. robustness in the geometry), cf. also Sogn:Takacs:2018 ().

In the present paper, we illustrate our findings with numerical experiments. All presented numerical experiments are available in the G+Smo library gismoweb ().

This paper is organized as follows. We give the model problem and a conforming discretization in Section 2. Then, in Section 3, we discuss why a non-conforming discretization might be of interest. Moreover, we propose a discontinuous Galerkin approach that fits our needs. We proceed to multigrid solvers: In Section 4, we discuss Gauss-Seidel smoothers and their performance. Motivated by that section, we introduce a subspace corrected mass smoother in Section 5 and finally a hybrid smoother in Section 6. In Section 7, we conclude and give some outlook. The Appending finally contains the proofs of the theorems stated in the paper.

## 2 Model problem and standard Galerkin discretization

Let with be an open and simply connected Lipschitz domain. Most of the numerical experiments are done for the two-dimensional domains shown in Figure 1.

The first domain is an L-shaped domain consisting of three quadrilaterals. Here, the geometry function is just the identity or a translation. On the coarsest grid level , each patch consists only of one element, i.e., the local basis functions are polynomials only. The second domain is the Yeti footprint which consists of the 21 patches depicted in Figure 1. Here, the grid on the coarsest grid level is as follows. The five patches at the bottom consist of two elements each, which are constructed by subdividing the patches on their longer sides. The remaining patches consist only of one element each. The grid levels are obtained by uniform refinement.

Consider the following Poisson model problem with Dirichlet boundary conditions. Find such that

 (1)

where and are given functions. Here and in what follows, , and are the standard Lebesgue and Sobolev spaces.

The experiments are performed for the choice and ; note that is the exact solution of the problem.

After homogenization (, ), the problem reads in variational form as follows. Find such that

 (∇u,∇v)L2(Ω)=(f,v)L2(Ω)for all v∈H10(Ω). (2)

The computational domain is a standard multi-patch domain. Thus, we assume that is composed of non-overlapping patches :

 ¯¯¯¯Ω=K⋃k=1¯¯¯¯¯¯Ωk with Ωk∩Ωl=∅ for k≠l, (3)

where each patch is represented by a sufficiently smooth bijective geometry function

 Gk:ˆΩ:=(0,1)d→Ωk:=Gk(ˆΩ)⊂Rd (4)

which can be continuously extended to , the closure of . Moreover, we assume that the mesh introduced by the patches satisfies the following condition.

###### Assumption 1

For any , the intersection is either (a) empty, (b) one common vertex, (c) the closure of one common edge, or – for – (d) the closure of one common face.

For each of the patches, we assume to have a hierarchy of grids with levels obtained by uniform refinement, which we denote by

 Vk,ℓ:={v∈L2(Ωk):v∘Gk∈d⨂δ=1Sp,hℓ}=span {φ(i)k,ℓ}Nk,ℓi=1, (5)

where is the space of tensor-product splines of degree , smoothness (or, equivalently, ) and grid size on the parameter domain . Note that the grid can be non-uniform and both the spline degree and the grid size can depend on the spatial direction and of the patch number; for simplicity, we do not write down this dependence explicitly. Note that the grid needs to be quasi-uniform, i.e., there needs to be a constant such that all knot spans on grid level are bounded from below by . The functions are assumed to form a (standard) B-spline or NURBS basis of .

To be able to set up a conforming discretization, we need to assume that the function spaces are fully matching on the interfaces, cf. (Takacs:2018, , Assumption 2.4). For tensor-product B-spline bases, the following assumption characterizes fully matching discretizations.

###### Assumption 2

On each interface between two patches, the geometry functions, the knot vector in tangential direction, and the spline degree in tangential direction agree.

Assuming a fully matching discretization, we define the conforming discretization space by

 Vcℓ:={v∈V:v|Ωk∈Vk,ℓ for k=1,…,K}. (6)

A basis for this space is visualized in Figure 2 (left), where all basis functions are represented by their Greville point. The support of the basis functions with Greville point in the interior of a patch is completely contained in that patch. The basis functions with Greville points on the interfaces are combinations of the matching patch-local basis functions. Their support extends to the vertices if and only if the Greville point is located on the vertex.

The conforming discretization of the model problem is obtained using the standard Galerkin principle: Find such that

 (uℓ,vℓ)H1(Ω)=(f,vℓ)L2(Ω)% for all vℓ∈Vcℓ. (7)

Using the abovementioned basis for the space , we obtain a standard matrix-vector problem: Find such that

 Aℓu––ℓ=f––ℓ, (8)

where is the stiffness matrix, the vector is the representation of with respect to the chosen basis and the load vector is obtained by testing the function with the basis functions.

## 3 Symmetric interior penalty discontinuous Galerkin (SIPG) discretization

Following LMMT:2015 (); LT:2015 (); Takacs:2019 (), we use a conforming isogeometric discretization for each patch and couple the patches using discontinuous Galerkin. We assume that the domain is again subdivided into patches such that (3), (4) and Assumption 1 are satisfied. We assume again to have patch-local spaces as in (5), which are combined in a non-conforming (i.e., discontinuous) way, i.e., we just define

 Vnℓ:={v∈L2(Ω):v|Ωk∈Vk,ℓ % for k=1,…,K and v|∂Ω=0}. (9)

This allows us to drop Assumption 2. Note that we strongly enforce the Dirichlet boundary conditions in our example; alternatively, one could use the SIPG method also to enforce the Dirichlet boundary conditions.

Since we have a discontinuous function space, we can visualize the degrees of freedom by tearing apart the individual patches, cf. Figure 3. Here, neither the Greville points nor the basis functions need to agree on the interfaces; the support of each basis function is contained in one single patch.

Since , it is not feasible to use the standard Galerkin principle for discretization. Thus, we couple the patches using the Symmetric Interior Penalty discontinuous Galerkin (SIPG) method. First, we define

 N:={(k,j):k

to be the set of interface-indices. For each interface with , we define the following symbols.

• n is the outer normal vector of . (Thus, is the outer normal vector of .)

• is the jump operator: .

• is the averaging operator: .

Now, we can formulate the SIPG discretization as follows: Find such that

 (uℓ,vℓ)Aℓ=(f,vℓ)L2(Ω)for all% vℓ∈Vnℓ, (10)

where we define

 (u,v)Aℓ :=(u,v)Qℓ−(u,v)Bℓ−(v,u)Bℓ, (u,v)Qℓ:=(u,v)Kℓ+σp2hL(u,v)Jℓ, (11) (u,v)Kℓ :=K∑k=1(∇u,∇v)L2(Ωk), (u,v)Jℓ:=∑(k,l)∈N(\llbracketu\rrbracket,\llbracketv\rrbracket)L2(Ik,l), (u,v)Bℓ :=∑(k,l)∈N(\llbracketu\rrbracket,{∇v}⋅{n})L2(Ik,l).

There is some independent of the grid size, the spline degree and the number of patches such that for all , the bilinear form is coercive, cf. (Takacs:2019, , Theorems 8 and 9). Thus, for , the Theorem of Lax Milgram states that the problem (10) has exactly one solution. The combination of Ceá’s Lemma and a naive approximation error estimate yields a discretization error estimate of the form

 |u−uL|2QL≤cp2h2L|u|2H2(Ω),

cf. Takacs:2019 (). By doing a more careful analysis, we obtain estimates of the form

 |u−uL|2QL≤c(logp)4h2L|u|2H2(Ω),

see (Takacs:2019, , eq. (15)). This significantly decreases the influence of the spline degree.

Note that the penalization term has the form

 σp2hL,

i.e., it depends on the grid size on the finest grid . This follows the ideas from Gopalakrishnan:Kanschat:2003 (). The idea behind that is that

 (uℓ,vℓ)Aℓ=(uℓ,vℓ)Aℓ+1and(uℓ,vℓ)Qℓ=(uℓ,vℓ)Qℓ+1 (12)

hold, i.e., we obtain a multigrid solver with conforming coarse-grid correction. This means that – on the coarse grid levels – the discretization is over penalized by a factor of , i.e.,

 σp2hL˜Σℓ:==2L−ℓσp2hℓΣℓ:=,

where is the canonical parameter and is the chosen one. We will see that this does not cause any problems for the examples we consider; following Gopalakrishnan:Kanschat:2003 (), convergence theory only holds if the number of smoothing steps is sufficiently increased for the coarser grid levels, cf. Remark 5.

Using a basis for the space , we obtain a standard matrix-vector problem: Find such that

 Aℓu––ℓ=f––ℓ. (13)

## 4 Multigrid solvers with Gauss-Seidel smoothers

In this and the following sections, we discuss several possible choices of multigrid smoothers, illustrate their convergence behavior with numerical experiments, and comment on the convergence theory.

We consider conforming discretizations and non-conforming discretizations which are set up as discussed in the last two sections. As we have nested spaces in all cases, the matrix is always the canonical embedding from into and the restriction matrix is its transpose. The chosen method is presented as pseudo-code as Algorithm 1, where we choose for the V-cycle method or for the W-cycle method.

In the finite element world, Gauss-Seidel smoothers are known to be very efficient smoothers; thus, as a first attempt, we consider such a smoother. One forward Gauss-Seidel sweep can be represented by

 u––ℓ←u––ℓ+L−1ℓ(f––ℓ−Aℓu––ℓ),

where is a lower-triangular matrix containing the coefficients of the stiffness matrix , i.e., it is given by

 (Lℓ)i,j={(Aℓ)i,j if i≥j0 if i

To be able to use our multigrid solver as preconditioner for a conjugate gradient solver, the post-smoothing procedure uses the transposed matrix , which represents a backward-Gauss-Seidel sweep.

All tables show the number of iterations required until the stopping criterion

 ∥ALu––L−f––L∥ℓ2∥f––L∥ℓ2≤ϵ:=10−8

is satisfied.

As usual, the convergence behavior of the overall solver can be improved if the multigrid method is not just used directly as a solver, but as a preconditioner within a preconditioned conjugate gradient (PCG) solver. Thus, we present results for both possibilities; in the following sections we will restrict ourselves to the more efficient PCG solver. Since the V-cycle and the W-cycle methods yield comparable iteration counts, we present the results for the more efficient V-cycle only. The number of smoothing steps is chosen as in all cases.

The multigrid solver was implemented in C++ based on the G+Smo library gismoweb (). The tables shown in the remainder of this section are obtained with the following command line instructions, where the values and are substituted accordingly.

> git clone https://github.com/gismo/gismo.git
> cd gismo
> make
> cd build/bin
> ./multiGrid_example -g domain2d/ldomain.xml -r  -p
-s gs -i d                                         // Table 1 (a)
> ./multiGrid_example -g domain2d/yeti_mp2.xml -r  -p
-s gs -i d                                         // Table 2 (a)

The results for the PCG experiments, presented in Tables 1 (b) and 2 (b), are obtained by replacing the option -i d by the option -i cg.

In Table 1, we observe that the multigrid solver is certainly robust in the grid size. While this approach is very efficient for small numbers of spline degrees, we observe that the convergence rates deteriorate significantly if the spline degree is increased. On the right side of the table, one can see the iteration counts for a preconditioned conjugate gradient method where one V-cycle of the mentioned multigrid method is used as a preconditioner. We observe that the iteration counts are significantly smaller than the iteration counts obtained by directly applying the multigrid solver. However, we simultaneously observe that we do not observe any qualitative improvement. In Table 2, we give the iteration counts for the Yeti footprint. We observe that – despite the fact that the geometry function is now non-trivial – the iteration counts are very similar to those of the L-shaped domain.

When one turns to the non-conforming discretizations, it immediately turns out that the multigrid solver utilizing the Gauss-Seidel smoother does not converge well at all.

One can show using standard arguments that the multigrid method converges with rates that are independent of the grid size and of the number of patches. The convergence analysis (for the conforming case) employs estimates that increase exponentially in the spline degree, cf. Gahalaut:Kraus:Tomar:2013 (). The numerical experiments show that this is not only a matter of the proof.

Since we did not obtain convincing results, we are interested in more advanced smoothers that work well also for SIPG discretizations and which do not deteriorate if the spline degree is increased.

## 5 Multigrid with subspace corrected mass smoother

In this section, we employ the subspace corrected mass smoother as introduced in Hofreither:Takacs:2017 (). That smoother requires that the spline space has a tensor-product structure; in our examples, we have such a structure on each patch but not on the overall domain. The extension of that smoother to conforming discretizations has been discussed in Takacs:2018 () based on a domain-decomposition approach. The key idea was to decompose all degrees of freedom on a per-piece bases. Pieces are the patch-interiors and the interface pieces. In two dimensions, the interface pieces are the edges and the vertices of each of the patches. In three dimensions, the interface pieces are the faces, the edges and the vertices of each of the patches. The decomposition of the degrees of freedom is depicted in Figure 4 (left).

The piece-local smoothers, which we denote by , are defined as follows. For being a patch-interior, we choose to be the subspace corrected mass smoother as proposed in Hofreither:Takacs:2017 (). We choose the scaling parameter (which was called  in Hofreither:Takacs:2017 ()) to be for some suitable chosen parameter . If is an interface piece, we choose

 Lℓ,T:=P⊤TAℓPT,

where the matrix represents the embedding of the piece in the whole space. The symbol refers to be the application of a direct solver. Applying a direct solver on the interfaces is feasible since the interfaces have much smaller numbers of degrees of freedom than the interiors of the patches. The overall smoother is just an additive composition of the piece-local smoothers, i.e., we choose

 L−1ℓ:=τ∑TPℓ,TL−1ℓ,TP⊤ℓ,T, (14)

where the sum is taken over all pieces . Here, is some damping parameter to be chosen.

The convergence theory from Takacs:2018 () can be summarized by the following theorem.

###### Theorem 3

Assume that is such that full elliptic regularity is satisfied (cf. (Takacs:2018, , Assumption 3.1)). Consider the conforming discretization and a multigrid solver with the smoother (14). There are constants , and which are independent of , , and (but may depend particularly on the geometry functions and the maximum number of neighbors of a patch) such that for

 τ∈(0,τ∗),δ∈(0,δ∗)andνℓ>ν∗ℓ:=pτ∗τδ∗δθ, (15)

the W-cycle multigrid method converges with a convergence rate .

Note that the terms and imply that the convergence degrades if too small values of those parameters are chosen. Thus, it is of interest to choose these parameters in a rather optimal way.

Note that the convergence theorem requires full elliptic regularity (cf. (Takacs:2018, , Assumption 3.1)). Thus, it is applicable to the Yeti footprint but it is not directly applicable to the L-shaped domain. Convergence results for the case with full elliptic regularity often carry over in practice to cases where that regularity assumption does not hold. The same behavior can be observed for the numerical experiments we have considered.

Observe that the convergence theorem suggests that the number of smoothing steps should increase with . As already outlined in Takacs:2018 (), this seems to be too pessimistic since the numerical experiments have shown that in all cases yields good convergence rates.

The numerical experiments are again applied within the same setup as in the last section. We set up a V-cycle multigrid method with 1+1 smoothing steps of the proposed smoother (on all grid levels). The damping parameter  is chosen as indicated with the option --MG.Damping and the scaling parameter  is chosen as indicated with the option --MG.Scaling below. The tables for the conforming case shown in this section are obtained with the following code, where the values and are substituted accordingly:

> ./multiGrid_example -g domain2d/ldomain.xml -r  -p        -s scms --MG.Damping 1 --MG.Scaling .12 -i d        // Table 3 (a) > ./multiGrid_example -g domain2d/yeti_mp2.xml -r  -p        -s scms --MG.Damping .25 --MG.Scaling .2 -i cg     // Table 5 (a)

The results for the PCG experiments presented Table 3 (b) are obtained by replacing the option -i d by the option -i cg.

All numerical experiments show that the proposed method is robust both in the grid size and the spline degree. However, when comparing the results for the Yeti footprint from Table 5 (a) with the corresponding results for the L-shaped domain from Table 3 (b), we see that the multigrid solver suffers from distorted geometry functions.

The numbers for the Yeti footprint seem not to be completely robust in the grid size. Since we have given a convergence theory, we know that the convergence numbers are bounded uniformly. Thus, the observed behavior is pre-asymptotic. The reason for this is that on coarser grid levels, the geometry is not resolved exactly. Let be the original stiffness matrix and be the simplified stiffness matrix obtained by neglecting the geometry function. Then, we have

 κ(ˆA−1ℓAℓ)=supvℓ∈Vℓ|vℓ|2H1(Ω)∑Kk=1|vℓ∘Gk|2H1(ˆΩ)supvℓ∈Vℓ∑Kk=1|vℓ∘Gk|2H1(ˆΩ)|vℓ|2H1(Ω),

which yields

 κ(ˆA−10A0) ≤⋯≤κ(ˆA−1L−1AL−1)≤κ(ˆA−1LAL) (16) ≤supv∈V|v|2H1(Ω)∑Kk=1|v∘Gk|2H1(ˆΩ)supv∈V∑Kk=1|v∘Gk|2H1(ˆΩ)|v|2H1(Ω),

which can be finally bounded by a constant times some power of the quantity . Of none of these relations is satisfied by equality. In such a case, the iteration counts are likely to increase if the grid gets refined. For more on this topic, see (Sogn:2018, , Section 7.4).

As a next step, we turn towards the non-conforming discretizations. Here, each degree of freedom is assigned to exactly one patch. So, it would be tempting to set up a patch-wise splitting of the degrees of freedom. Unfortunately, numerical experiments have shown that this approach does not work well. So, we follow the approach from Takacs:2018 () also in the non-conforming case and split the degrees of freedom again into pieces . This means that we avoid breaking the coupling which was enforced by the penalty term. So, the degrees of freedom belonging to one edge (face, vertex) are considered to be one piece, even if the degrees of freedom belong to different patches, see Figure 4 (right).

For this choice, we can give the following convergence theorem.

###### Theorem 4

Assume that is such that full elliptic regularity holds (cf. (Takacs:2019, , Assumption 4)) and assume that the geometry functions (but not necessarily the discretizations) agree on the interfaces (cf. (Takacs:2019, , Assumption 2)). Consider the SIPG discretization and a multigrid solver with the smoother (14). There are constants , and which are independent of , , and (but may depend particularly on the geometry functions and the maximum number of neighbors of a patch) such that for

 τ∈(0,τ∗),δ∈(0,δ∗)andνℓ>ν∗ℓ:=2L−ℓ(L−ℓ)2p(logp)4τ∗τδ∗δθ, (17)

the W-cycle multigrid method converges with a convergence rate .

We give the proof of this theorem in the Appendix; the proof is based on the error estimates from Takacs:2019 ().

One might observe that the number of smoothing steps required by this convergence theorem increases like . This follows the approach suggested in Gopalakrishnan:Kanschat:2003 () and is related to the chosen over-penalization discussed in Section 3. {remark} Note that the number of degrees of freedom on the coarser grid levels is smaller by a factor of . So, also using these additional smoothing steps, the overall complexity of the multigrid solver is still linear in the number of unknowns on the finest grid level if (a) or (b) the V-cycle is considered. If we consider and the W-cycle, the choice (17) yields that the computational complexity grows like , where is the number of unknowns on the finest grid level. In the numerical experiments, we did not observe that increasing the number of smoothing steps has been required. Analogously to the conforming case, also the stated dependence on is too pessimistic; thus, we again choose on all grid levels.

Now, we provide numerical experiments for the SIPG discretization. Theoretically, we could just use exactly the discretization that has been chosen for the conforming case. This, however, yields a (particularly uninteresting) special case since Assumption 2 holds. In this special case, we have and the SIPG formulation converges to the conforming discretization for . Instead, we are interested in a discretization such that Assumption 2 does not hold: We modify the setup of the spaces. For one third of the patches, we use the original spline space . For one third of the patches, we use the spline space . For the last third of the patches, we use the spline space . This particular setting is obtained with the command line option --NonMatching. In this way, we obtain a setup where a conforming discretization is not possible.

The tables for the non-conforming case shown in this section are obtained with the following code, where the values and are substituted accordingly:

> ./multiGrid_example -g domain2d/ldomain.xml -r  -p  --DG       --NonMatching -s scms --MG.Damping .9 --MG.Scaling .12       -i d                                            // Table 4 (a) > ./multiGrid_example -g domain2d/yeti_mp2.xml -r  -p  --DG       --NonMatching -s scms --MG.Damping .25 --MG.Scaling .2       -i cg                                           // Table 5 (b)

The results for the PCG experiments presented in Table 4 (b) are obtained by replacing the option -i d by the option -i cg.

The PCG discretizations are presented in Tables 4 (b) and 5 (b); we again obtain robustness in the grid size and the spline degree. Here, for the Yeti footprint, we have again iteration counts that are slightly increasing with the grid size; again, this observation can be explained by the fact that finer grids allow to resolve the geometry function better, cf. (16).

In principle, the method works also if the multigrid solver is applied directly, cf. Table 4 (a). Here, we suffer from numerical instabilities which are amplified with an increasing number of levels. One can avoid these instabilities, e.g., by increasing the number of pre- and post smoothing steps. However, using the multigrid method as a preconditioner within a PCG solver is obviously the more efficient approach.

## 6 Multigrid with hybrid smoother

We have observed that a multigrid method with the subspace corrected mass smoother is robust in the grid size and the spline degree and works well for both conforming and discontinuous Galerkin discretizations. We have also observed that this approach suffers from non-simple geometry functions since it is based on the close connection between the stiffness matrix and the simplified stiffness matrix . The results for the Gauss-Seidel smoother are different: the multigrid solver works badly both for large spline degrees and for discontinuous Galerkin discretizations. However, by comparing Table 1 with Table 2, we observe that the method behaves quite robust in the geometry function.

Since the behavior of the two smothers is somewhat orthogonal, we can hope for a good method if we combine them. Our idea is to use one forward Gauss-Seidel sweep followed by the subspace corrected mass smoother for pre-smoothing and the subspace corrected mass smoother followed by one backward Gauss-Seidel sweep for post-smoothing. The overall method is presented as Algorithm 2.

The convergence analysis from Section 5 can be easily carried over to the hybrid smoother. The iteration matrix for the (V or W cycle) multigrid method with the hybrid smoother is given by

 ˜Wℓ:=(I−(LGSℓ)−⊤Aℓ)Wℓ(I−(LGSℓ)−1Aℓ),

where is the iteration matrix of the (V or W cycle, respectively) multigrid method with the subspace corrected mass smoother. Since the Gauss-Seidel iteration is stable in the energy norm, we obtain

 ∥˜Wℓ∥Aℓ≤∥I−(LGSℓ)−⊤Aℓ∥Aℓ∥Wℓ∥Aℓ∥I−(LGSℓ)−1Aℓ∥Aℓ≤∥Wℓ∥Aℓ.

So, we have using the results from the last section the convergence of the W-cycle multigrid method with hybrid smoother. Thus, we obtain as follows.

###### Corollary 5

Consider the multigrid solver with the hybrid smoother. Under the assumptions of Theorem 3 or 4, respectively, the W-cycle multigrid method converges with a convergence rate .

The tables for the experiments with the hybrid smoother are obtained with the following code, where the values and are substituted accordingly:

> ./multiGrid_example -g domain2d/yeti_mp2.xml -r  -p        -s hyb --MG.Damping .25 --MG.Scaling .1 -i d       // Table 6 (a) > ./multiGrid_example -g domain2d/yeti_mp2.xml -r  -p  --DG       --NonMatching -s hyb --MG.Damping .25 --MG.Scaling .1       -i d                                               // Table 7 (a)

The results for the PCG experiments, presented in Tables 6 (b) and 7 (b), are obtained by replacing the option -i d by the option -i cg.

For both cases, we obtain that the iteration counts are quite robust in the grid size (even if the maximum number of iterations is not reached for the coarser grid levels). We observe that the number of iterations increases with the spline degree. This is indeed due to the fact that for small values of , the Gauss-Seidel smoother yields very fast convergence and that convergence behavior is carried over to the hybrid smoother. For larger spline degrees, the hybrid smoother’s convergence behavior degrades mildly; this is due to the fact that the Gauss-Seidel smoother is not completely capable to capture all effects perfectly. Still, keeping in mind that the condition number of the stiffness matrix grows exponentially with the spline degree, the observed behavior is still very satisfactory.

Compared to applying the subspace corrected mass smoother only, the hybrid smoother pays of in most cases. Certainly, applying the hybrid smoother with means basically that pre- and post-smoothing steps are applied. Since the Gauss-Seidel smoother is slightly cheaper than the subspace corrected mass smoother, the costs for one such cycle are smaller than the costs of two multigrid cycles with the subspace corrected mass smoother only.

Besides the two-dimensional examples considered so far, the proposed methods can be directly extended to three dimensional problems (even if the details of the convergence theory have not been worked out for these cases). We consider the two domains depicted in Figure 5: the Fichera corner and a variant of that domain with non-trivial geometry function, which we call twisted Fichera corner.

The tables for the three dimensional domains are obtained with the following code, where the values and are substituted accordingly:

> ./multiGrid_example -g domain2d/fichera.xml -r  -p        -s hyb --MG.Scaling .12 --MG.Damping 1 -i cg        // Table 8 (a) > ./multiGrid_example -g domain2d/twisted_fichera.xml -r  -p        -s hyb --MG.Scaling .12 --MG.Damping .25 -i cg     // Table 9 (a)

The results for the DG experiments, presented in the Tables 8 (b) and 9 (b), are obtained by adding the command line options --DG --NonMatching.

Similar to the results for the Yeti footprint, Tables 8 and 9 again show small iteration counts. For the twisted Fichera corner, we observe that the number of iterations increases mildly when the grid gets refined; this is again to be explained by the better resolution of the geometry function. Moreover, we observe a very mild dependence on the spline degree.

## 7 Conclusions and outlook

We have presented robust multigrid solvers for multi-patch IgA with conforming and non-conforming discretizations. We have given convergence results that exactly state the robustness of the solvers in the grid size. Concerning the dependence on the spline degree, the statements seem to be too pessimistic since the solvers have been completely or (at least) rather robust in practice.

We have addressed another issue, which causes problems for all solvers that use the tensor-product structure on the parameter domain: the dependence on the geometry function. We have proposed a hybrid smoother between our subspace corrected mass smoother and the Gauss-Seidel smoother which seems to reduce the effect on the geometry function. Finding approaches to better incorporate the geometry function into the smoother itself seems to be an interesting topic for further research.

## Appendix

In the appendix, we give a proof of Theorem 4 and of some auxiliary results.

Every constant used within the appendix is assumed to be independent of the grid size, the grid level, the spline degree and the number of patches, but it may depend on the geometry function (cf. (Takacs:2019, , Assumption 3)), the number of neighbors of a patch (cf. (Takacs:2018, , Assumption 2.3)), the constant in the elliptic regularity assumption (cf. (Takacs:2018, , Assumption 3.1)) and the quasi-uniformity of the grid, i.e., the ratio between the largest and the smallest knot span of the knot vectors on one level. We write if and only if there is a constant such that and we write if and only if and .

First, we show the following lemma, which is basically a trace inequality.

###### Lemma 6

Let be the space of tensor-product splines of degree on a quasi-uniform grid with size on the parameter domain . Then, the estimate

holds for all and all .

{proof}

For , let be the eigenfunctions of , i.e., such that

 (ψν,i,ψν,j)L2(0,1)=δi,jand (ψ′ν,i,ψ′ν,j)L2(0,1)+θ2(ψν,i,ψν,j)L2(0,1)=λν,iδi,j,

where is the Kronecker delta and are the corresponding eigenvalues. Using coercivity of and a standard inverse estimate, cf. (Schwab:1998, , Corollary 3.94), we obtain

 θ2≤λν,1andλν,Nν≲p4h−2+θ2.

We define level sets

 Iν,m:={i:μm−1:=2m−1θ2≤λν,i<μm:=2mθ2}

for , where

 M:=1+maxν∈{1,2}⌊log2(θ−2λν,Nν)⌋≲log(1+p4h−2θ−2)

is the number of level sets. Note that by construction every eigenvalue belongs to exactly one level set. Every function can be represented as

 u(x,y) =N1∑i=1N2∑j=1ui,jψ1,i(x)ψ2,j(y)=M∑m=1M∑n=1∑i∈I1,m∑j∈I2,nui,jψ1,i(x)ψ2,j(y)wm,n(x,y):=.

A standard trace estimate, cf. (Takacs:2018, , Lemma 4.4), yields

 |wm,n(0)|2 ≲∥wm,n∥L2({0}×(0,1))∥wm,n∥H1({0}×(0,1)) (18) ≲∥wm,n∥1/20,0,1∥wm,n∥1/21,0,1∥wm,n∥1/20,1,1∥wm,n∥1/21,1,1

where

 ∥w∥2a,b,η:=∥∥∥∂a+b∂xa∂xbw∥∥∥2L2(ˆΩ)+η2(a+b)∥w∥2L2(ˆΩ)

for . Since and since all eigenvalues are in or , respectively, we obtain

 ∥w∥2a,b,1≤∥w∥2a,b,θ≂μamμbn∥wm,n∥2L2(ˆΩ).

Using (18), we obtain further

 |wm,n(0)|2 ≲μmμn∥wm,n∥2L2(ˆΩ)≂∥wm,n∥1,0,θ∥wm,n∥0,1,θ ≤(|wm,n|2H1(ˆΩ)+θ2∥wm,n∥2L2(ˆΩ)).

Finally, the Cauchy-Schwarz inequality and orthogonality of the basis functions (both in and ) yield

 |u(0)|2 ≲M2M∑m=1M∑n=1|wm,n(0)|2 ≲M2M∑m=1M∑n=1(|wm,n|H1(ˆΩ)+θ2∥wm,n∥L2(ˆΩ))2 =M2(|u|2H1(ˆΩ)+θ2∥u∥2L2(ˆΩ)),

which finishes the proof. ∎

Now, we give bounds on the smoother which allow to show the smoothing property.

###### Lemma 7

Provided the assumptions of Theorem 4, the estimate

 Aℓ≤Lℓ≲p(logp)2(L−ℓ)22L−ℓτ∗τδ∗δ˜Lℓ

holds, where and is the standard mass matrix.

{proof}

The proof of this Lemma requires the notation from Takacs:2018 (), i.e., we denote the set of all patch-interiors by , the set of all edges by and the set of all vertices by .

Observe that we have

 Aℓ≂Qℓ (19)

for all , where is defined by (11) and (12). For , this statement directly follows from (Takacs:2019, , Theorem 8). Since (Takacs:2019, , Theorem 8) also holds in cases of over-penalization, we can apply that theorem also to the case and obtain (19) also in that case.

As first step, we bound from below. (The following arguments are analogous to (Takacs:2018, , Lemma 4.3).) The triangle inequality yields

 Aℓ≲∑T∈K∪E∪VPℓ,T(P⊤ℓ,TAℓP⊤ℓ,T)P⊤ℓ,T. (20)

For , (Hofreither:Takacs:2017, , Lemma 8) and (19) yield . For , we have by definition . Thus, we obtain from (20)

 Aℓ≲∑T∈K∪E∪VPℓ,TLℓ,TP⊤ℓ,T

and for all with small enough further

 Aℓ≤τ−1∑T∈K∪E∪VPℓ,TLℓ,TP⊤ℓ,T=Lℓ,

which shows the first part of the desired inequality.

Now, we bound from above. We use the decomposition

 Qℓ=Kℓ+σp2hLJℓ,

cf. (11). Using (19), we obtain

 Lℓ,T=P⊤ℓ,TAℓPℓ,T≂P⊤ℓ,TQℓPℓ,T=P⊤ℓ,TKℓPℓ,T˜Kℓ,T:=+σp2hLP⊤ℓ,TJℓPℓ,T˜Jℓ,T:=

for all and, therefore,

 Lℓ≂˜Kℓ+σp2hL˜Jℓ, (21)

where

 ˜Kℓ:=τ−1∑T∈KPℓ,TLℓ,TP⊤ℓ,T+τ−1∑T∈E∪VPℓ,TKℓ,TP⊤ℓ,T

and

 ˜Jℓ:=τ−1∑T∈E∪VPℓ,TJℓ,TP⊤ℓ,T.

Completely analogous to (Takacs:2018, , Lemma 4.7), we obtain

 ˜Kℓ≲pτ∗τδ∗δ(Kℓ+h−2ℓMℓ)≤pτ∗τδ∗δ(Qℓ+h−2ℓMℓ)≤p2L−ℓτ∗τδ∗δ˜Lℓ. (22)

So, it remains to bound from above. Since the restriction of to any patch-interior vanishes, the same arguments as in the proof of (Takacs:2018, , Lemma 4.7) and the triangle inequality yield

 ∑T∈E</