- Email: [email protected]
Overlapping Schwarz preconditioners for isogeometric collocation methods
Accepted Manuscript Overlapping Schwarz preconditioners for isogeometric collocation methods L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi PII: DOI: Reference:
S0045-7825(14)00152-2 http://dx.doi.org/10.1016/j.cma.2014.05.007 CMA 10241
To appear in:
Comput. Methods Appl. Mech. Engrg.
Received date: 22 March 2014 Revised date: 12 May 2014 Accepted date: 15 May 2014 Please cite this article as: L.B. da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, Overlapping Schwarz preconditioners for isogeometric collocation methods, Comput. Methods Appl. Mech. Engrg. (2014), http://dx.doi.org/10.1016/j.cma.2014.05.007 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Overlapping Schwarz preconditioners for Isogeometric collocation methods L. Beir˜ao da Veigaa , D. Chob , L. F. Pavarinoa , S. Scacchia a
b
Dipartimento di Matematica, Universit` a di Milano, Via Saldini 50, 20133 Milano, Italy Department of Mathematics, Dongguk University, Pil-dong 3-ga, Jung-gu, Seoul, 100-715, South Korea.
Abstract Isogeometric collocation methods are very recent and promising numerical schemes that preserve the advantages of isogeometric analysis but often exhibit better performances than their Galerkin counterparts. In the present paper, an additive overlapping Schwarz method for isogeometric collocation discretizations is introduced and studied. The resulting preconditioner, accelerated by GMRES, is shown to be scalable with respect to the number of subdomains and very robust with respect to the isogeometric discretization parameters such as the mesh size and polynomial degree, as well as with respect to the presence of discontinuous elliptic coefficients and domain deformations. 1. Introduction Isogeometric Analysis (IGA) is an innovative numerical methodology for the solution of partial differential equations (PDEs), introduced by Hughes et al. in [20], that potentially allows for a direct connection with CAD, thus providing a much easier and exact representation of the computational domain in a wide range of applications. Indeed, in IGA the same basis functions representing the CAD geometry (for instance NURBS or T-splines) are also used for the discretization of the unknown fields, following an isoparametric paradigm. Since its introduction, IGA techniques have had a lot of success and have been studied and applied in diverse fields; see the monograph [14] and the first mathematical paper [5] for a general introduction to this field. As described in the contributions above, since its first introduction IGA has been mainly following a Galerkin approach, i.e. adopting a discretization of the variational form of the problem. More recently it has been observed that, due to high continuity of the discrete space and the presence of non-polynomial functions, numerical integration in IGA (necessary in order to build the stiffness matrix) can be particularly costly when compared to the number of degrees of freedom [21, 4]. Moreover, the large advantage in terms of accuracy granted by the high continuity can be somewhat offset by the denser structure of the associated matrices, see for instance [12, 13]. Collocation methods have recently been introduced in IGA as a class of schemes that allow to overcome the above limitations, see the papers [1, 2]. Indeed, building the collocation stiffness matrix requires only one evaluation per basis function and, moreover, the resulting linear system turns out to be much sparser than in Galerkin methods. A cost comparison between collocation and Galerkin IGA schemes can be found in [29]. Collocation IGA methods have also been used very successfully in beam analysis, see [9, 3], and for the Cahn-Hilliard problem [25]. Clearly all this advantages come at a price, that is the lack of a theoretical background in more than one dimension and the fact that these collocation schemes are not guaranteed to preserve some energy Email addresses: [email protected] (L. Beir˜ ao da Veiga), [email protected] (D. Cho), [email protected] (L. F. Pavarino), [email protected] (S. Scacchi) Preprint submitted to Computer Methods in Applied Mechanics and Engineering
May 22, 2014
principles. Note moreover that, even for self-adjoint problems, the stiffness collocation matrixes are, in general, not symmetric. As in the Galerkin case, the collocation stiffness matrix has a condition number that grows very rapidly when the mesh is refined and/or the polynomial degree is raised. As a consequence, preconditioners become a very important tool, since otherwise the high performance of collocation methods in terms of cost-to-accuracy can be lost because of (1) the accuracy degeneration due to machine precision algebra combined with a bad stability constant and (2) the large number of iterations needed for iterative solvers. Preconditioners and domain decomposition methods [31, 30] have a recent but successful history in Galerkin schemes for IGA, see for instance [6, 7, 8, 10, 26, 11, 23, 22], while there is no result in the literature for preconditioners or domain decomposition methods for IGA collocation schemes. In the present contribution, we develop and analyze an overlapping additive Schwarz (OAS) domain decomposition method for IGA collocation schemes. We focus on the following model elliptic problem: −∇ · (ρ∇u) = f in Ω, u = 0 on ∂Ω, (1) where the scalar field 0 < ρmin ≤ ρ(x) ≤ ρmax for all x ∈ Ω and Ω ⊂ Rd is a bounded and connected CAD domain. The proposed scheme is based on a decomposition of the basis functions index set into overlapping subsets, associated to “rectangular” overlapping subdomains in parameter space. Once these sets are introduced, a one-level preconditioner can be built in an additive way by solving local collocation problems on the overlapping subdomains. Since one-level preconditioners cannot be expected to be scalable with the number of subdomains, we propose also a two-level preconditioner that employs an additional coarse problem. A total of three different coarse problems are proposed. The first two employ a nested coarse space, one obtained by a projection of the fine matrix on the coarse space, and the other obtained by enforcing the equations at a new set of coarse collocation points. The third coarse problem follows the same approach as the first one above, but employs a non-nested coarse space. After presenting our isogeometric OAS preconditioner, we test its performance through a series of numerical tests in two and three dimensions. Since the matrices of collocation schemes are not symmetric, we accelerate our preconditioner with GMRES [32] in order to solve iteratively the resulting collocation problem. We first show the low number of iterations required to solve the OAS preconditioned system (as opposed to the unpreconditioned system), studying how the iteration counts depend on the mesh size and number of subdomains. In particular, from the numerical results one can appreciate the scalability of the two-level OAS method. We also study the method behavior in terms of the polynomial degree p and regularity index k. The iteration counts are very robust with respect to these parameters, even for minimal overlap, a property that is not shared by minimal overlap OAS in Galerkin IGA. We also test the scheme behavior in presence of domain deformation and discontinuous elliptic coefficients, with very satisfactory results. We remark that, although the numerical tests in the current contribution have been implemented on a single-core machine, the proposed OAS preconditioner is also immediately suitable for a parallel implementation on multi-core machines. The present paper is organized as follows. In Section 2 we give a brief presentation of NURBS and isogeometric analysis. In Section 3 we describe the proposed overlapping domain decomposition preconditioner. In Section 4 we present the results of extensive numerical tests in two and three dimensions.
2
2. NURBS-based isogeometric analysis Non-Uniform Rational B-splines are a standard tool for describing and modeling curves and surfaces in computer aided design and computer graphics (see Piegl and Tiller [27] and Rogers [28] for an extensive description of these functions and their properties). In this work, we use NURBS as an analysis tool, as proposed by Hughes et al. [20]. In this section, we present a short description of B-splines, NURBS, the basics of isogeometric analysis and an introduction to the proposed discretization method. In the rest of the paper, we will adopt the following compact notation. Given two real numbers a, b we write a . b, when a ≤ Cb for a generic constant C independent of the knot vectors (defined below), and we write a ≈ b when a . b and b . a. 2.1. B-splines B-splines in the plane are piecewise polynomial curves composed of linear combinations of Bspline basis functions. A knot vector is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve {ξ1 = 0, ..., ξn+p+1 = 1},
(2)
where p is the polynomial degree of the B-spline and n is the number of basis functions (and control points) necessary to describe it. The open interval (ξ1 , ξn+p+1 ) is called a patch. A knot vector is said to be uniform if its knots are uniformly-spaced and non-uniform otherwise. The maximum allowed knot multiplicity is p + 1; a knot vector is said to be open if its first and last knots have multiplicity p + 1. In what follows, we always employ open knot vectors. Basis functions formed from open knot vectors are interpolatory at the ends of the parametric interval Ib := (0, 1) but are not, in general, interpolatory at interior knots. Given a knot vector, univariate B-spline basis functions are defined recursively starting with p = 0 (piecewise constants) ( 1 if ξi ≤ ξ < ξi+1 Ni0 (ξ) = (3) 0 otherwise. For p > 1 : ξ − ξi N p−1 (ξ) + ξi+p+1 − ξ N p−1 (ξ) if ξ ≤ ξ < ξ i i+p+1 p ξi+p+1 − ξi+1 i+1 Ni (ξ) = ξi+p − ξi i 0 otherwise,
(4)
where, in (4), we adopt the convention 0/0 = 0. Thus, the general basis function Nip has support Θi := supp(Nip ) = (ξi , ξi+p+1 ),
i = 1, 2, .., n.
The functions Nip are a partition of unity, see e.g. [14]. In Figure 1 we present an example consisting of n = 9 cubic basis functions generated from the open knot vector ξ = {0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1} If internal knots are not repeated, B-spline basis functions are C p−1 -continuous. If a knot has multiplicity α, the basis is C p−α -continuous at that knot. In particular, when a knot has multiplicity p, the basis is C 0 and interpolates the control point at that location. In the following, we will assume that the maximum knot multiplicity is p so that all considered functions will be (at least) globally continuous. We define the spline space Sbh = span{Nip (ξ), i = 1, . . . , n}. 3
(5)
By means of tensor products, a multi-dimensional B-spline region can be constructed. We discuss here the case of a two-dimensional region, the higher-dimensional case being analogous. b := (0, 1) × (0, 1) be the two-dimensional parametric space. Consider the knot vectors {ξ1 = Let Ω 0, ..., ξn+p+1 = 1} and {η1 = 0, ..., ηm+q+1 = 1}, and an n × m net of control points Ci,j . One dimensional basis functions Nip and Mjq (with i = 1, ..., n and j = 1, ..., m) of degree p and q, b is then defined by respectively, are defined from the knot vectors. The bivariate spline basis on Ω the tensor product construction p,q Bi,j (ξ, η) = Nip (ξ) Mjq (η).
Observe that the two knot vectors {ξ1 = 0, ..., ξn+p+1 = 1} and {η1 = 0, ..., ηm+q+1 = 1} generate a mesh of rectangular elements in the parametric space in a natural way. Analogous to (5), we define p,q Sbh = span{Bi,j (ξ, η), i = 1, . . . , n, j = 1, . . . , m}
(6)
2.2. NURBS maps and spaces
In general, a rational B-spline in Rd is the projection onto d-dimensional physical space of a polynomial B-spline defined in (d+1)-dimensional homogeneous coordinate space. For a complete discussion, see the book by Farin [19] and references therein. In this way, a great variety of geometrical entities can be constructed and, in particular, all conic sections in physical space can be obtained exactly. To obtain a NURBS curve in R2 , we start by introducing the NURBS basis functions of degree p N p (ξ)ωi N p (ξ)ωi Rip (ξ) = Pn i p = i , (7) w(ξ) ˆı=1 Nˆı (ξ)ωˆı P where the denominator w(ξ) = ˆnı=1 Nˆıp (ξ)ωˆı ∈ Sbh is called the weight function. The NURBS curve is then defined by C(ξ) =
n X
Rip (ξ)Ci ,
(8)
i=1
where Ci ∈ R2 are control points. Analogously to B-splines, NURBS basis functions on the two-dimensional parametric space b = (0, 1) × (0, 1) are defined as Ω p,q Ri,j (ξ, η)
p,q p,q Bi,j (ξ, η)ωi,j Bi,j (ξ, η)ωi,j = Pn Pm = , p,q w(ξ, η) ˆı=1 ˆ=1 Bˆı,ˆ (ξ, η)ωˆı,ˆ
(9)
where ωi,j = (Cωi,j )3 and the denominator is the weight function denoted also by w(ξ, η). Observe that the continuity and support of NURBS basis functions are the same as for B-splines. NURBS spaces are the span of the basis functions (9). NURBS regions are defined in terms of the basis functions (9). In particular a single-patch domain Ω is a NURBS region associated with the n × m net of control points Ci,j , and we introduce b → Ω given by the geometrical map F : Ω F(ξ, η) =
n X m X i=1 j=1
4
p,q Ri,j (ξ, η)Ci,j .
(10)
Following the isoparametric approach, the space of NURBS scalar fields on the domain Ω is defined, component by component as the span of the push-forward of the basis functions (9) p,q Nh := span{Ri,j ◦ F−1 , with i = 1, . . . , n; j = 1, . . . , m}.
(11)
As observed for instance in [5], in order to obtain spaces with homogeneous Dirichlet boundary conditions it is sufficient to eliminate the first and last function in each coordinate. In the present paper, in order to keep simple homogeneous Dirichlet boundary conditions but make use of a general structure, we do not exploit such a possibility. Therefore the isogeometric space used in the discretization will be Vb = Sbh and V = Nh living in parameter space and in physical space, respectively and the boundary conditions will be enforced in the discrete problem. The interested reader may find more details on isogeometric analysis as well as many interesting applications in a number of recently published papers mentioned in the introduction. The threedimensional case is analogous and not discussed here. 2.3. Isogeometric collocation In the present section we briefly review the isogeometric collocation method for the present problem (1), introduced in [1], see also [2, 29]. We focus mainly on the bi-dimensional case, the one and three dimensional cases being defined analogously. For simplicity of exposition, we consider directly the case where the collocation points are chosen as the Greville abscissae associated to the knot vectors, see [16]. The choice of collocation points is crucial for the stability and good behavior of the discrete problem; another possible (but less common) choice is given by the Demko abscissae [17]. Let ξ i , i = 1, ..., n, be the Greville abscissae related to the knot vector {ξ1 , ..., ξn+p+1 }: . ξ i = (ξi+1 + ξi+2 + ... + ξi+p )/p .
(12)
Analogously, let η j , j = 1, ..., m, be the Greville abscissae related to the knot vector {η1 , ..., ηm+q+1 }. It is easy to see that ξ 1 = η 1 = 0, ξ n = η m = 1, while all the remaining points are in (0, 1). We define the collocation points τij ∈ Ω by the tensor product structure τij = F(b τij ) ,
τbij = (ξ i , η j ) ∈
b , Ω
for i = 1, ..., n, j = 1, ..., m. Then, the isogeometric collocation problem with Greville abscissae reads: find uh ∈ V such that ( − ∇ · (ρ∇uh )(τij ) = f (τij ) i = 2, ..., n − 1, j = 2, ..., m − 1 , (13) uh (τij ) = 0 (i, j) ∈ ({1, n} × {1, ..., m}) ∪ ({1, ..., n} × {1, m}) . The well-posedness and convergence to the exact solution for problem (13) has been studied, only for the one-dimensional case, in [1]. The theory for the isogeometric collocation method in two and three dimensions is still an open issue; nevertheless, a large array of numerical tests in the literature seem to indicate that the method is stable and convergent in all practical cases. Note that we consider homogeneous Dirichlet boundary conditions only for simplicity; different kind of boundary conditions can be clearly treated in the same fashion.
5
Remark 2.1. In order to ease the presentation of the collocation method we here assume that the knot repetitions are such that the ensuing discrete spline and NURBS spaces are at least C 2 , so that the equations in (13) are well defined. In other words, the polynomial degree p ≥ 3 and the number of repetitions of any knot (excluding the boundary ones) never exceeds p − 2. Should this assumption fail, the method can be easily modified as detailed in [2]; we show an example in Section 4.5. 3. The overlapping Schwarz preconditioner In this section, we construct an isogeometric overlapping Schwarz preconditioner for the isogeometric collocation method (13). 3.1. Subdomains and subspace decomposition We describe first the subdomains and subspace decompositions in 1D and then extend them by tensor products to 2D and 3D. The decomposition is first built for the underlying space of spline functions in parameter space, and then easily extended to the NURBS space in the physical domain. We select from the full set of knots {ξ1 = 0, ..., ξn+p+1 = 1} a subset {ξik , k = 1, . . . , N +1} of (non-repeated) interface knots with ξi1 = 0, ξiN+1 = 1. This subset of interface knots defines a decomposition of the closure of the reference interval Ib = [0, 1] =
[
k=1,..,N
Ibk , with Ibk = (ξik , ξik+1 ),
that we assume to have a similar characteristic diameter H ≈ Hk = diam(Ibk ). The interface knots are thus given by ξik for k = 2, .., N . For each of the interface knots ξik we choose an index 2 ≤ sk ≤ N − 1 (strictly increasing in k) that satisfies sk < ik < sk + p + 1, so that the support of the basis function Nspk intersects both Ibk−1 and Ibk . Note that at least one such sk exists; if it is not unique, any choice can be taken. N
3 i
N
1
0
3 i
1
0
0.2
0.4
ξ
0.6
0.8
0
1
(a) r = 0
0
0.2
0.4
ξ
0.6
0.8
1
(b) r = 1
Figure 1: Cubic basis functions formed from ξ = {0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1}. For various r, Vb1 is the span of basis functions drawn with dash-dot and solid lines and Vb2 is the span of basis functions drawn with solid and b In particular, the basis functions in common dashed lines, in two subdomains Ib1 = (0, 1/2) and Ib2 = (1/2, 1) of I. are those drawn with a solid line. Small rectangles at the bottom denote the Greville abscissae related to the knot vector ξ.
6
We then define an overlapping decomposition of Ib in the following way. Let r ∈ N be an integer (called the overlap index) counting the basis functions shared by adjacent subdomains. We define the index sets k = 1, 2, .., N, Θk = {j ∈ N : sk − r ≤ j ≤ sk+1 + r}
with the exception that 1 ≤ j ≤ s2 + r for the index set Θ1 and sN − r ≤ j ≤ n for the index set ΘN and the local spaces Vbk = span{Njp (ξ) : j ∈ Θk }
k = 1, 2, .., N,
(14)
with the analogous exception for Vb1 , VbN . An example is shown in Figure 1. These subspaces form an overlapping decomposition of the spline space Vb ; indeed, 2r + 1 represents the number of basis functions in common (in the univariate case) among “adjacent” local subspaces. For r = 0 we have the minimal overlap consisting of just one common basis function between subspaces. We now define the extended subdomains Ibk′ by Ibk′ =
[
Njp ∈Vbk
supp(Njp ) = (ξsk −r , ξsk+1 +r+p+1 ),
(15)
′ . Let us then introduce a (open) coarse knot vector with the analogous exception for Ib1′ and IbN 0 ξ 0 = {ξ10 = 0, ..., ξN = 1} corresponding to a coarse mesh determined by the subdomains Ibk , c +p+1 i.e. ξ 0 = {ξ1 , ξ2 , . . . , ξp , ξi1 , ξi2 , ξi3 , . . . , ξiN−1 , ξiN , ξiN +1 , ξiN +2 . . . , ξiN +p+1 }, (16)
such that the distance between adjacent distinct knots is of order H, ξ1 = · · · = ξp = ξi1 = 0 and ξiN+1 = ξiN +2 = · · · = ξiN +p+1 = 1 and the associated coarse spline space Vb0 := SbH = span{Ni0,p (ξ), i = 1, ..., Nc }
has the same degree p as Sbh and is thus a subspace of Sbh . Note that, in order to have a smaller coarse space Vb0 , in the definition above we repeat only once all the internal knots (see also Remark 3.1). In more than one dimensions, we just proceed by tensor products. For example, in 2D we define subdomains and overlapping subdomains by Ibk = (ξik , ξik+1 ), b ′ = Ib′ × Ib′ , Ω kl
k
l
Ibl = (ηjl , ηjl+1 ),
b kl = Ibk × Ibl Ω
1 ≤ k ≤ N, 1 ≤ l ≤ M.
N M Moreover, we take the indices {sk }N k=2 associated to the {ξik }k=2 and the analogous indices {sl }l=2 M associated to the {ηjl }l=2 . The local index sets are defined by
Θkl = {(i, j) ∈ N2 : sk − r ≤ i ≤ sk+1 + r, sl − r ≤ j ≤ sl+1 + r} for k = 1, 2, .., N and l = 1, 2, .., M . The local and coarse subspaces now are p,q Vbkl = span{Bi,j (ξ, η) : (i, j) ∈ Θkl }, ◦ p,q
◦ p,q
Vb0 = span{B i,j : B i,j (ξ, η) := Ni0,p (ξ)Mj0,q (η), i = 1, ..., Nc , j = 1, ..., Mc }. 7
The decomposition for the NURBS space V in the physical domain is the trivial extension of the one above. Therefore the local subspaces and the coarse space are, up to the usual modifications for the boundary subdomains, p,q Vkl = span{Ri,j ◦ F−1 : (i, j) ∈ Θkl }, ◦ p,q
◦ p,q
V0 = span{Ri,j ◦ F−1 := (B i,j /w) ◦ F−1 , i = 1, ..., Nc , j = 1, ..., Mc }, where we recall that w is the weight function, see (9). The subdomains in physical space are defined as the image of the subdomains in parameter space with respect to the mapping F b kl ), Ωkl = F(Ω
b ′ ). Ω′kl = F(Ω kl
Remark 3.1. Note that less regular (but still nested) coarse spaces can be considered by including more repetitions of the interior knots. The only requirement for nestedness is that each interior knot is repeated at most as many times as the corresponding knot in the vector defining the fine space. Although we here focus on the more regular case (since it leads to the smaller nested coarse space) the cases with knot repetitions are handled exactly in the same way in what follows, see (18) and (19). 3.2. Matrix form of the preconditioner We are now able to present the preconditioner Schwarz operator, focusing as usual on the bidimensional case. Let the (local and coarse) restriction matrixes Rkl : V → Vkl and R0 : V → V0 be T : V → V , k = 1, .., N , l = 1, .., M defined as the transpose of the natural embedding matrixes Rkl kl T R0 : V0 → V . Let Akl , for all k = 1, .., N , l = 1, .., M , be the square matrix associated to the local collocation problems: find ukl h ∈ Vkl such that ( − ∇ · (ρ∇ukl (i, j) ∈ Θikl , h )(τij ) = f (τij ) ukl h (τij ) = 0
(i, j) ∈ Θ∂kl ,
where the internal and boundary index sets are Θikl = (i, j) ∈ Θkl : (i, j) ∈ / ({1, n} × {1, ..., m}) ∪ ({1, ..., n} × {1, m}) Θ∂kl = Θkl \Θikl .
Note that the above matrix exactly corresponds to T Akl = Rkl A Rkl
(17)
where A is the global collocation matrix associated to the original problem (13). Choice of coarse problem. We propose three choices for the coarse problem. The first two make use of a nested coarse space while the third one makes use of a non-nested coarse space. • First nested choice (A0 ). We simply define, in analogy to (17), the coarse matrix by A0 = R0 A R0T ,
(18)
where R0T is the matrix form of the interpolation operator from the coarse space V0 to the fine space V . Note that the above matrix corresponds to a collocation problem in which each equation in the system is obtained as a weighted sum of the collocation equations of the differential equation on more than one point. 8
• Second nested choice (A′0 ). We consider the collocation matrix ensuing from the coarse basis functions and associated Greville abscissae. We thus define A′0 as the square matrix associated to the coarse collocation problem: find u0h ∈ V0 such that ( − ∇ · (ρ∇u0h )(τij0 ) = f (τij0 ) i = 2, ..., Nc − 1, j = 2, ..., Mc − 1 , u0h (τij0 ) = 0
(i, j) ∈ ({1, Nc } × {1, ..., Mc }) ∪ ({1, ..., Nc } × {1, Mc }), (19) 0 where τij are the Greville abscissae associated to the coarse knot vector (16) by the usual rule (12). In what follows, A′0 denotes the global collocation matrix associated to the coarse problem (19).
• Non-nested choice (A′′0 ). The coarse matrix is defined as in (18), but in this case the coarse space is constituted by the bi-linear or tri-linear piecewise polynomials on the coarse mesh. Since V0 is not a subspace of V , R0T is the matrix form of the L2 projection from V0 to V . Finally, the matrix form of the proposed overlapping Schwarz operator is TOAS = BOAS A, where BOAS is the additive Schwarz preconditioner BOAS = R0T A−1 0 R0 +
N X M X
T Rkl A−1 kl Rkl .
(20)
k=1 l=1
The coarse matrix A0 is to be replaced by the matrix A′0 or A′′0 according to the choice of coarse problem. The three dimensional case is handled analogously. e kl assoMore general multiplicative and hybrid versions as well as projection-like operators T ciated with inexact local solvers could also be used, see [30, 31], but will not be considered here. The use of the BOAS preconditioner (20) for the iterative solution of the discrete system Au = f can also be regarded as replacing it with the preconditioned system TOAS u = g, where g = BOAS f , which can be accelerated by a Krylov subspace method. Remark 3.2. Note that, as is standard in collocation methods, the matrix A is in general not symmetric, even if the original problem is self-adjoint. The same holds for the preconditioned operator; therefore GMRES (instead of CG - conjugate gradient) iterations will need to be used in the solution process. On the other hand, note that the assembly of the local and coarse matrixes Akl , A0 is much cheaper than in Galerkin IGA, exactly as is the case for the global matrix A (see [29]). Moreover such matrixes will be sparser than their Galerkin counterparts. The numerical results presented in the next Section indicate that the convergence rate of our collocation OAS solver is comparable to those of the analogous Galerkin OAS solver studied in our previous work [6], i.e. it is independent of the number of subdomains N (scalability), optimal in the ratio H/h, robust with respect to the IGA polynomial degree p and regularity index k, as well as with respect to the presence of discontinuous elliptic coefficients and domain deformations. Unlike in the IGA Galerkin case, a mathematical analysis of these good convergence rate properties is still missing in the IGA collocation case, as well as a mathematical analysis of the more basic approximation properties of IGA collocation. 9
h 1/8 1/16 1/32 1/64
p=3 29.19 28.56 82.10 327.21
h 1/8 1/16 1/32 1/64
p=3 21.67 84.15 334.15 1334.23
IGA Galerkin p=4 p=5 p=6 240.03 2148.25 20818.97 222.55 1812.58 15132.64 215.00 1700.63 13634.81 381.73 1688.11 13493.15 IGA collocation p=4 p=5 p=6 36.51 54.94 177.74 138.42 205.88 287.38 546.52 810.23 1126.91 2179.25 3228.47 4486.68
p=7 222889.58 129178.94 110155.73 108472.87 p=7 695.15 589.37 1497.23 5955.86
Table 1: Condition numbers computed with respect to the 2-norm of the IGA Galerkin and IGA collocation stiffness matrices on the unit square domain, as a function of the mesh size h and of the spline polynomial degree p. Spline regularity is always maximal, thus k = p − 1.
4. Numerical results In the present section, we test the convergence properties of the isogeometric overlapping Schwarz preconditioner (20) for 2D and 3D model elliptic problems (1) on both parametric (reference square or cube) and physical domains, with zero right-hand side, Dirichlet boundary conditions ex sin(y) in 2D or mixed boundary conditions in 3D. The problem is discretized with isogeometric NURBS spaces with associated mesh size h, polynomial degree p, regularity k, using as a starting basis the Matlab isogeometric library GeoPDEs [15]. The domain is decomposed into N overlapping subdomains of characteristic size H and overlap index r, as described in Section 3.1. The discrete problems are solved by the GMRES method with the isogeometric Schwarz preconditioner (20), with a zero initial guess and as stopping criterion a 10−6 reduction of the relative residual. In the following tests, we study how the convergence rate of the OAS preconditioner depends on h, N, p, k, on the jumps of the elliptic coefficients and on domain deformations. We will mostly consider as collocation points the tensor-product Greville abscissae (12), but in some tests we will also consider the Demko abscissae [17] to show that our OAS preconditioner performs equally well for other choices of collocation points. 4.1. 2D tests: condition number of the collocation IGA stiffness matrix In this test, we study the behavior of the condition number of the IGA collocation stiffness matrix with respect to the mesh size h and the spline polynomial degree p. The domain is the unit square and the regularity of the splines is maximal, thus k = p − 1. The condition number is computed with respect to the 2-norm. The results reported in Table 1 show that the condition number of the IGA collocation stiffness matrix increases as h−2 , as that of the IGA Galerkin stiffness matrix for p = 3. We observe that, for coarse meshes and high p, the condition number of the IGA Galerkin stiffness matrix initially decreases. This unexpected behavior has already been observed in [24]. Up to p = 7, the growth with respect to p of the IGA collocation stiffness matrix seems to be much weaker than that of the IGA Galerkin stiffness matrix, which has an exponential growth, see [24].
10
N 2×2 4×4 8×8 16 × 16 32 × 32
* 12
1/h = 8 A0 A′′0 13 13 -
A′0 13
p = 3, k = 2, r = 0 square domain 1/h = 16 1/h = 32 A0 A′′0 A′0 * A0 A′′0 A′0 * 14 15 15 22 16 19 17 29 16 19 21 29 18 21 25 40 36 17 20 52 53 69 -
* 17 21
1/h = 64 A0 A′′0 A′0 21 25 19 22 25 26 18 22 61 16 19 66 -
* 39 55 73 97 167
1/h = 128 A0 A′′0 A′0 26 34 23 28 34 25 23 25 62 18 21 89 16 19 83
Table 2: 1-level (*) and 2-level OAS preconditioners with nested (A0 ), non-nested (A′′0 ) and second choice (A′0 ) coarse spaces, on the unit square domain: GMRES iteration counts as a function of the number of subdomains N and mesh size 1/h. Fixed p = 3, k = 2, r = 0.
35
21 A’’0, p=4, k=3
20
A’’0, p=3, k=2
19
30 0
GMRES ITERATIONS
GMRES ITERATIONS
A’’ , p=3, k=2 18
A0, p=4, k=3
17 16
A0, p=3, k=2
15 14 13
A’’ , p=4, k=3 0
25
A , p=3, k=2 0
A , p=4, k=3 0
20
15
12 11 10 0
200
400 600 800 NUMBER OF SUBDOMAINS N
1000
10 0
1200
10
20
30
40 50 RATIO H/h
60
70
80
Figure 2: GMRES iteration counts for increasing number of subdomains and fixed H/h = 4 (left panel) and for increasing ratio H/h and fixed N = 2 × 2 (right panel). 2-level Overlapping Schwarz preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse spaces, minimal overlap r = 0, (p = 3, k = 2) and (p = 4, k = 3) spline spaces on the unit square domain.
N 2×2 4×4 8×8 16 × 16 32 × 32
1/h = 8 * A0 A′′0 13 13 12 -
p = 3, k = 2, 1/h = 16 * A0 A′′0 17 14 15 21 16 18 -
r = 0 square domain 1/h = 32 1/h = 64 ′′ * A0 A0 * A0 A′′0 22 16 19 29 20 25 30 18 21 40 22 25 36 17 20 52 19 22 69 16 20 -
1/h = 128 * A0 A′′0 38 27 35 54 28 33 73 24 26 100 18 22 161 16 19
Table 3: 1-level (*) and 2-level OAS preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse spaces, on the unit square domain: GMRES iteration counts as a function of the number of subdomains N and mesh size 1/h. Fixed p = 3, k = 2, r = 0. Demko abscissae are considered as collocation points.
11
N 2×2 4×4 8×8 16 × 16 32 × 32
p = 3, k = 2, r = 0 quarter of ring domain 1/h = 8 1/h = 16 1/h = 32 1/h = 64 * A0 A′′0 * A0 A′′0 * A0 A′′0 * A0 A′′0 16 15 14 21 18 17 28 22 22 37 28 28 28 21 21 39 26 24 55 34 31 53 26 22 76 33 26 104 30 23 -
1/h = 128 * A0 A′′0 50 36 38 73 45 39 109 42 33 171 37 26 245 31 22
Table 4: 1-level (*) and 2-level OAS preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse space in the quarter of ring domain: GMRES iteration counts as a function of the number of subdomains N and mesh size 1/h. Fixed p = 3, k = 2, r = 0.
4.2. 2D tests: OAS scalability in N and optimality in H/h The GMRES iteration counts of the 1-level and 2-level OAS preconditioners on the unit square are reported in Table 2 for fixed p = 3, k = 2, r = 0. Additional results for p = 4, k = 3, r = 0 are plotted in Figure 2. Moreover, Table 3 shows that numerical results obtained by using as collocation points the Demko abscissae [17] are very similar to those in Table 2, indicating that our OAS preconditioners are robust with respect to the choice of collocation points. In what follows, the Greville abscissae are chosen as collocation points. The GMRES iteration counts of the 1-level and 2-level OAS preconditioners on the ring-shaped physical domain for fixed p = 3, k = 2, r = 0 are reported in Table 4. In all the three tables, the iteration counts are reported as function of the number of subdomains N and mesh size 1/h. The results show that the OAS(2) preconditioner is scalable with both nested and non-nested coarse spaces, since moving along the diagonal of the tables the OAS(2) iteration counts are bounded from above by a constant independent of N . This is shown clearly in the left panel of Figure 2. As expected, the coarse problem is essential for scalability, since the results corresponding to the 1-level preconditioner (i.e. without coarse problem) show iteration counts growing with N along the diagonals of the tables. The results also show that the iteration counts depend less than linearly on the ration H/h (right panel of Figure 2 or moving along rows in the tables). These results for varying N and H/h indicate that the proposed collocation OAS preconditioners performs as well as the analogous OAS preconditioners we have previously investigated for Galerkin IGA discretizations (see [6, Sec. 5.1]). 4.3. 2D test: complex spectrum of the OAS preconditioned operator In order to investigate the good convergence rate observed in the previous numerical tests, we report in Figure 3 the spectrum of the 2D OAS preconditioned operator TOAS = BOAS A in the complex plane, considering three cases with increasing number of subdomains N = 16, 64, 144 (first, second, third column, respectively). The increasing ill-conditioning of the unpreconditioned matrix is shown by the spectra of A (top panels), since for increasing N the eigenvalues seem to approach zero and grow unbounded in modulus (for N = 16 and 64 the spectrum of A turns out to be real). On the other hand, the spectrum of the 2-level preconditioner (bottom panels) is contained in a complex region [α, 5] × [−1, 1] with α > 0.5. This fact guarantees the existence of a sequence of complex polynomials with value 1 at the origin and decreasing values on the OAS spectrum, i.e. a fast GMRES convergence, see e.g. [32]. The spectrum of the 1-level preconditioned operator (middle panels) is contained in a similar region but where α seems to approach 0 when N increases, making it impossible to find a sequence of complex polynomials with the same property as before, thus confirming the non-scalability of OAS(1). 12
4
4
2
10
0
2
10
0
0
50 100 150 200 250 min real=9.81e+00 max real=1.39e+03 min im=0.00e+00 max im=0.00e+00
10
300
eig B−1A, OAS(1)
x 10
1 0
eig B
−1 −2
0
0
200 400 600 800 1000 min real=9.88e+00 max real=5.72e+03 min im=0.00e+00 max im=0.00e+00
1200
0
−4
−1
A, OAS(1)
−15
2
0.2 0.1 0 −0.1 −0.2
0.5 1 1.5 2 2.5 3 3.5 4 min real=1.62e−01 max real=4.00e+00 min im=−7.22e−03 max im=7.22e−03
4
x 10
2 0 −2 −4
4.5
0
0.5 1 1.5 2 2.5 3 3.5 4 min real=4.42e−02 max real=4.00e+00 min im=−4.40e−03 max im=4.40e−03
4
eig B
−0.005 −0.01
0
0.5 1 1.5 2 2.5 3 3.5 4 4.5 min real=7.15e−01 max real=4.66e+00 min im=−7.22e−03 max im=7.22e−03
5
x 10
0
4.5
0.5 1 1.5 2 2.5 3 3.5 4 4.5 min real=6.30e−01 max real=4.75e+00 min im=−3.83e−03 max im=3.83e−03
5
0
0 −3
0
−5
0.5 1 1.5 2 2.5 3 3.5 4 min real=1.99e−02 max real=4.00e+00 min im=−3.83e−03 max im=3.83e−03
x 10
−2 −4
4.5
eig B−1A, OAS(2)
eig B−1A, OAS(2)
0
−1
A, OAS(2)
0.005
5
14000
2
−3
0.01
2000 4000 6000 8000 10000 12000 min real=9.90e+00 max real=1.31e+04 min im=−2.05e−01 max im=2.05e−01 −4
eig B−1A, OAS(1)
10
eig A
10
eig A
eig A
10
0.5 1 1.5 2 2.5 3 3.5 4 4.5 min real=6.54e−01 max real=4.73e+00 min im=−4.40e−03 max im=4.40e−03
5
5
x 10
0
−5
0
Figure 3: Spectrum of the collocation unpreconditioned (first row), 1-level (second row) and 2-level (third row) OAS preconditioned operators. 2D quarter-ring domain, N = 16 (first column), 64 (second column), 144 (third column), H/h = 4, p = 3, k = 2.
N k=2 k=3 k=4 k=5 k=6
h = 1/64, N = 4 × 4(H/h = 16), r p=3 p=4 p=5 NP A0 A′′0 NP A0 A′′0 NP A0 215 22 25 307 20 26 414 17 316 23 23 399 22 393 20 -
= 0 square domain p=6 A′′0 NP A0 A′′0 26 467 26 30 28 433 17 26 25 444 21 30 513 20 22 -
NP 655 648 668 473 570
p=7 A0 30 18 18 18 19
A′′0 32 31 26 26 25
Table 5: Unpreconditioned (NP) GMRES and 2-level OAS preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse space, on the unit square domain: GMRES iteration counts as a function of the spline polynomial degree p and the regularity k. Fixed 1/h = 64, N = 4 × 4, H/h = 16, symmetric minimal overlap (r = 0).
13
ρ
NP
OAS(1)
1e-4 1e-2 1e+0 1e+2 1e+4
> 104 1753 496 > 104 > 104
37 48 70 58 62
ρ=1
α
ρ = 10
OAS(2) A0 A′′0 39 39 32 42 30 42 50 42 87 45
Table 6: Central jump test on the unit square domain (left panel). OAS robustness with respect to jump discontinuities in the elliptic coefficient ρ (right table): Iteration counts of the unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with both nested A0 and non-nested A′′0 coarse spaces) preconditioned GMRES. h = 1/64, N = 8 × 8, H/h = 8, overlap r = 0, p = 3, k = 2.
4.4. 2D tests: OAS dependence on p and k The aim of this test is to compare the GMRES iteration counts of the unpreconditioned case and 2-level OAS preconditioners (with nested or non-nested coarse spaces) with respect to the polynomial degree p and spline regularity k. The domain considered is the unit square discretized by a mesh of size h = 1/64, subdivided into N = 4 × 4 subdomains with symmetric minimal overlap r = 0. In this case, if k is odd, we take two common functions instead of one between adjacent univariate subspaces in order to preserve symmetry in terms of basis functions. Table 5 reports the iteration counts, varying p from 3 to 7 and k from 2 to 6. The iteration counts of the unpreconditioned case increase when p increases, while the iteration counts of the 2-level OAS preconditioner seem to have a very robust behavior with respect to k and p. In particular, we note that such behavior is markedly better than its counterpart in [6] for Galerkin IGA discretizations. 4.5. 2D tests: OAS robustness with respect to jumping coefficients In order to investigate the robustness of the OAS preconditioners with respect to jumps discontinuities of the elliptic coefficient ρ, we consider the classical “central jump” test in the unit square decomposed into 8 × 8 subdomains. The elliptic coefficient ρ varies by 8 orders of magnitude (from 10−4 to 104 ) in the 4 × 4 central subdomains, while it is 1 in the surrounding subdomains, see the left figure in Table 6. We fix h = 1/64, H/h = 8, p = 3, k = 2, r = 0. The regularity across the jump interface is reduced to C 0 . Following [2], we enforce the continuity of the weighted normal derivative at the collocation points on the interface layer. Table 6 reports the iteration counts for the unpreconditioned (NP) system, 1-level OAS (OAS(1)) and 2-level OAS (OAS(2)) preconditioned GMRES, with both nested (A0 ) and non-nested (A′′0 ) coarse spaces. The results show clearly the robustness of 1-level OAS and 2-level OAS with non-nested coarse space, since the GMRES iterations remain bounded when ρ jumps. On the contrary, the nested case seems to be more sensitive to high values of the jumps. 4.6. 2D tests: OAS robustness with respect to domain deformations We now study how the OAS preconditioners perform in presence of domain deformations. The four domains in Figure 4 are taken into account, deforming the rectangle [0, 0.5] × [0, 2] (A) into increasingly more curved boomerang-shaped domains (B, C, D). All the domains are discretized by the same mesh of size h = 1/64, subdivided into N = 4× 4 subdomains. The NURBS discretization space is (p = 3, k = 2). Both Greville and Demko collocation points are considered. The results show that both the 1-level and 2-level OAS preconditioners are considerably less sensitive to domain 14
2
2
2
2
1.5
1.5
1.5
1.5
1
1
1
1
0.5
0.5
0.5
0.5
0 0
0.5
(a)
0 0
0.5
1
0 0
(b)
0.5
1
1.5
0 0
0.5
(c)
1
1.5
2
2.5
(d)
Figure 4: Domains for the 2D boomerang test.
domain A B C D
Greville abscissae NP OAS(1) OAS(2) A0 A′′0 194 46 34 32 223 49 39 36 290 50 45 43 492 52 52 51
Demko abscissae NP OAS(1) OAS(2) A0 A′′0 177 45 32 31 214 50 37 36 302 50 45 44 485 52 53 52
Table 7: 2D boomerang test: OAS robustness with respect to domain deformations. Iteration counts of unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with nested A0 and non-nested A′′0 coarse spaces) preconditioned GMRES with Greville and Demko abscissae. h = 1/64, N = 4 × 4, H/h = 16, overlap r = 0, p = 3, k = 2.
deformations than the unpreconditioned solver. Table 7 shows that the GMRES iteration counts increase (going from A to D) by a factor 2.5 for the unpreconditioned solver, while only by a factor 1.5 for 2-level preconditioners. We note that the number of iterations for the 1-level preconditioner seems almost independent of the domain deformation. Nevertheless, since OAS(1) is non-scalable, we expect its iteration counts to increase for increasing subdomains N , while we expect the OAS(2) iteration counts to remain bounded. 4.7. 3D tests: OAS scalability in N The iteration counts of unpreconditioned (NP) GMRES, 1-level OAS (OAS(1)) and 2-level OAS with nested coarse space (OAS(2)) preconditioned GMRES are reported in Table 8 for the reference cubic domain, as a function of the number of subdomains N for fixed subdomain size H/h = 4 (i.e. both h and H are decreasing proportionally as in a scaled speedup test), p = 3, k = 2 and overlap r = 0. These 3D results confirm that 2-level OAS is scalable, since the condition number seems bounded from above by a constant independent of N . 4.8. 3D tests: OAS robustness with respect to domain deformations We finally study how the OAS preconditioners perform in presence of 3D domain deformations. The four domains in Figure 5 are taken into account, deforming the parallelepipedal domain [0, 0.5]× [0, 0.5] × [0, 2] (A) into increasingly more curved boomerang-shaped domains (B, C, D). All the 15
NP 24 33 44 52 64
N 2×2×2 3×3×3 4×4×4 5×5×5 6×6×6
OAS(1) 17 23 29 33 39
OAS(2), A0 17 22 25 26 26
Table 8: OAS scalability in N on the unit cube. Iteration counts of unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with nested A0 coarse space) preconditioned GMRES as a function of the number of subdomains N . Fixed H/h = 4, p = 3, k = 2.
2.5 2 1.5
1.5 1
1
1
0.5
0.5
0.5
0.5
0 2
0 2
0 2
0 2
1.5
1.5
1
1 0.5
0.5
0.5
0.5 0
0
0
(a)
1.5
1
1
(b)
2.5 1.5
1.5
1 0.5
0.5 0 0
0
(c)
2 1.5
1
1
0.5 0
0.5 0
(d)
Figure 5: Domains for the 3D boomerang test.
domains are discretized by the same mesh of size h = 1/12, subdivided into N = 3 × 3 × 3 subdomains. The NURBS discretization space has parameters p = 3, k = 2. Both Greville and Demko collocation points are considered. The results show that, as in the 2D case, both the 1-level and 2-level OAS preconditioners are less sensitive to domain deformations than the unpreconditioned solver. Table 9 shows that the GMRES iteration counts increase by a factor 3 (going from A to D) for the unpreconditioned solver, while only by a factor 1.3 for the 1-level and 2-level preconditioners. References [1] F. Auricchio, L. Beir˜ao da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli. Isogeometric Collocation Methods. Math. Mod. Meth. Appl. Sci., 20 (11), 2075–2107, 2010. [2] F. Auricchio, L. Beir˜ao da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli. Isogeometric collocation for elastostatics and explicit dynamics. Comp. Meth. Appl. Mech. Engrg., 249–252: 2–14, 2012. [3] F. Auricchio, L. Beir˜ao da Veiga, J. Kiendl, C. Lovadina, A. Reali. Locking-free isogeometric collocation methods for spatial Timoshenko rods. Comp. Meth. Appl. Mech. Engrg., 263: 113– 126, 2013.
16
domain A B C D
Greville abscissae NP OAS(1) OAS(2) A0 A′′0 57 32 30 39 68 33 33 33 99 35 37 37 179 38 41 40
Demko abscissae NP OAS(1) OAS(2) A0 A′′0 56 32 30 30 65 33 34 34 91 34 37 36 157 39 41 40
Table 9: 3D boomerang test: OAS robustness with respect to domain deformations. Iteration counts of unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with nested A0 and non-nested A′′0 coarse spaces) preconditioned GMRES with Greville and Demko abscissae. h = 1/12, N = 3 × 3 × 3, H/h = 4, overlap r = 0, p = 3, k = 2.
[4] F. Auricchio, F. Calabro, T.J.R. Hughes, A. Reali and G. Sangalli. A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis. Comput. Methods Appl. Mech. Engrg., 249-252: 15–27, 2012. [5] Y. Bazilevs, L. Beir˜ao da Veiga, J.A. Cottrell, T.J.R. Hughes, G. Sangalli. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. Math. Mod. Meth. Appl. Sci., 16, 1–60, 2006. [6] L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. Overlapping Schwarz methods for Isogeometric Analysis. SIAM J. Numer. Anal., 50:1394–1416, 2012. [7] L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. Isogeometric Schwarz preconditioners for linear elasticity systems. Comp. Meth. Appl. Mech. Engrg., 253:439–454, 2013. [8] L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. BDDC preconditioners for Isogeometric Analysis. Math. Mod. Meth. Appl. Sci., 23:1099–1142, 2013. [9] L. Beir˜ao da Veiga, C. Lovadina, A. Reali. Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods. Comp. Meth. Appl. Mech. Engrg., 241–244:38– 51, 2012. [10] L. Beir˜ao da Veiga, L.F. Pavarino, S. Scacchi, O.B. Widlund, S. Zampini. Isogeometric BDDC Preconditioners with Deluxe Scaling. To appear in SIAM J. Sci. Comput. [11] A. Buffa, H. Harbrecht, A. Kunoth, G. Sangalli. BPX-preconditioning for isogeometric analysis. Comp. Meth. Appl. Mech. Engrg., 265:63–70, 2013. [12] N. Collier, D. Pardo, L. Dalcin, M. Paszynsk, V.M. Calo. The cost of continuity: a study of the performance of isogeometric finite elements using direct solvers. Comp. Meth. Appl. Mech. Engrg., 213–216:353–361, 2012. [13] N. Collier, L. Dalcin, D. Pardo and V.M. Calo. The cost of continuity: Performance of iterative solvers on isogeometric finite elements. SIAM J. Sci. Comput., 35, 767–784, 2013. [14] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs. Isogeometric Analysis. Towards integration of CAD and FEA. Wiley, 2009. [15] C. De Falco, A. Reali, R. Vazquez. GeoPDEs: A research tool for Isogeometric Analysis of PDEs. Adv. Engrg. Softw., 42, 1020–1034, 2011. 2010. 17
[16] C. de Boor. A practical guide to splines. Springer, 2001. [17] S. Demko. On the existence of interpolation projectors onto spline spaces, J. of Approx. Theory, 43, 151–156, 1985. [18] M. Dryja, O. B. Widlund. Domain Decomposition Algorithms with Small Overlap, SIAM J. Sci.Comp., 15(3), 604–620, 1994. [19] G.E. Farin. NURBS curves and surfaces: from projective geometry to practical use. A.K. Peters, 1995. [20] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement. Comp. Meth. Appl. Mech. Engrg., 194, 4135– 4195, 2005. [21] T.J.R. Hughes, A. Reali, G. Sangalli. Efficient Quadrature for NURBS-based Isogeometric Analysis. Comp. Meth. Appl. Mech. Engrg., 199, 301–313, 2010. [22] K.P.S. Gahalaut, J.K. Kraus, S.K. Tomar. Multigrid methods for isogeometric discretization. Comp. Meth. Appl. Mech. Engrg., 253:413–425, 2013. [23] K.P.S. Gahalaut, S.K. Tomar, J.K. Kraus. Algebraic multilevel preconditioning in isogeometric analysis: Construction and numerical studies . Comp. Meth. Appl. Mech. Engrg., 266:40–56, 2013. [24] K.P.S. Gahalaut and S.K. Tomar. Condition number estimates for matrices arising in the isogeometric discretizations. Tech. Report RICAM, 23, 2012. [25] H. Gomez, A. Reali, G. Sangalli. Accurate, efficient, and (iso)geometrically flexible collocation methods for phase-field models. Accepted for publication on J. Comp. Phys. [26] S.K. Kleiss, C. Pechstein, B. J¨ uttler, S. Tomar. IETI-Isogeometric Tearing and Interconnecting. Comp. Meth. Appl. Mech. Engrg., 247–248:201–215, 2012. [27] L. Piegl, W. Tiller. The NURBS Book, 2nd Edition. Springer-Verlag, 1997. [28] D.F. Rogers. An Introduction to NURBS With Historical Perspective. Academic Press, 2001. [29] D. Schillinger, J. A. Evans, A Reali, M. A. Scott, T. J.R. Hughes. Isogeometric Collocation: cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations. Comp. Meth. Appl. Mech. Engrg., 267:170–232, 2013. [30] B. F. Smith, P. Bjørstad, W. D. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 1996. [31] A. Toselli, O. B. Widlund. Domain Decomposition Methods: Algorithms and Theory. Computational Mathematics, Vol. 34. Springer-Verlag, Berlin, 2004. [32] L. N. Trefethen, D. Bau. Numerical Linear Algebra. SIAM, 1996.
18
HIGHLIGHTS • We study isogeometric collocation methods which are very recent and promising numerical schemes. • We present an overlapping additive isogeometric collocation discretizations.
Schwarz
method
for
• We show that the overlapping additive Schwarz preconditioner is scalable and very robust.
S0045-7825(14)00152-2 http://dx.doi.org/10.1016/j.cma.2014.05.007 CMA 10241
To appear in:
Comput. Methods Appl. Mech. Engrg.
Received date: 22 March 2014 Revised date: 12 May 2014 Accepted date: 15 May 2014 Please cite this article as: L.B. da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, Overlapping Schwarz preconditioners for isogeometric collocation methods, Comput. Methods Appl. Mech. Engrg. (2014), http://dx.doi.org/10.1016/j.cma.2014.05.007 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Overlapping Schwarz preconditioners for Isogeometric collocation methods L. Beir˜ao da Veigaa , D. Chob , L. F. Pavarinoa , S. Scacchia a
b
Dipartimento di Matematica, Universit` a di Milano, Via Saldini 50, 20133 Milano, Italy Department of Mathematics, Dongguk University, Pil-dong 3-ga, Jung-gu, Seoul, 100-715, South Korea.
Abstract Isogeometric collocation methods are very recent and promising numerical schemes that preserve the advantages of isogeometric analysis but often exhibit better performances than their Galerkin counterparts. In the present paper, an additive overlapping Schwarz method for isogeometric collocation discretizations is introduced and studied. The resulting preconditioner, accelerated by GMRES, is shown to be scalable with respect to the number of subdomains and very robust with respect to the isogeometric discretization parameters such as the mesh size and polynomial degree, as well as with respect to the presence of discontinuous elliptic coefficients and domain deformations. 1. Introduction Isogeometric Analysis (IGA) is an innovative numerical methodology for the solution of partial differential equations (PDEs), introduced by Hughes et al. in [20], that potentially allows for a direct connection with CAD, thus providing a much easier and exact representation of the computational domain in a wide range of applications. Indeed, in IGA the same basis functions representing the CAD geometry (for instance NURBS or T-splines) are also used for the discretization of the unknown fields, following an isoparametric paradigm. Since its introduction, IGA techniques have had a lot of success and have been studied and applied in diverse fields; see the monograph [14] and the first mathematical paper [5] for a general introduction to this field. As described in the contributions above, since its first introduction IGA has been mainly following a Galerkin approach, i.e. adopting a discretization of the variational form of the problem. More recently it has been observed that, due to high continuity of the discrete space and the presence of non-polynomial functions, numerical integration in IGA (necessary in order to build the stiffness matrix) can be particularly costly when compared to the number of degrees of freedom [21, 4]. Moreover, the large advantage in terms of accuracy granted by the high continuity can be somewhat offset by the denser structure of the associated matrices, see for instance [12, 13]. Collocation methods have recently been introduced in IGA as a class of schemes that allow to overcome the above limitations, see the papers [1, 2]. Indeed, building the collocation stiffness matrix requires only one evaluation per basis function and, moreover, the resulting linear system turns out to be much sparser than in Galerkin methods. A cost comparison between collocation and Galerkin IGA schemes can be found in [29]. Collocation IGA methods have also been used very successfully in beam analysis, see [9, 3], and for the Cahn-Hilliard problem [25]. Clearly all this advantages come at a price, that is the lack of a theoretical background in more than one dimension and the fact that these collocation schemes are not guaranteed to preserve some energy Email addresses: [email protected] (L. Beir˜ ao da Veiga), [email protected] (D. Cho), [email protected] (L. F. Pavarino), [email protected] (S. Scacchi) Preprint submitted to Computer Methods in Applied Mechanics and Engineering
May 22, 2014
principles. Note moreover that, even for self-adjoint problems, the stiffness collocation matrixes are, in general, not symmetric. As in the Galerkin case, the collocation stiffness matrix has a condition number that grows very rapidly when the mesh is refined and/or the polynomial degree is raised. As a consequence, preconditioners become a very important tool, since otherwise the high performance of collocation methods in terms of cost-to-accuracy can be lost because of (1) the accuracy degeneration due to machine precision algebra combined with a bad stability constant and (2) the large number of iterations needed for iterative solvers. Preconditioners and domain decomposition methods [31, 30] have a recent but successful history in Galerkin schemes for IGA, see for instance [6, 7, 8, 10, 26, 11, 23, 22], while there is no result in the literature for preconditioners or domain decomposition methods for IGA collocation schemes. In the present contribution, we develop and analyze an overlapping additive Schwarz (OAS) domain decomposition method for IGA collocation schemes. We focus on the following model elliptic problem: −∇ · (ρ∇u) = f in Ω, u = 0 on ∂Ω, (1) where the scalar field 0 < ρmin ≤ ρ(x) ≤ ρmax for all x ∈ Ω and Ω ⊂ Rd is a bounded and connected CAD domain. The proposed scheme is based on a decomposition of the basis functions index set into overlapping subsets, associated to “rectangular” overlapping subdomains in parameter space. Once these sets are introduced, a one-level preconditioner can be built in an additive way by solving local collocation problems on the overlapping subdomains. Since one-level preconditioners cannot be expected to be scalable with the number of subdomains, we propose also a two-level preconditioner that employs an additional coarse problem. A total of three different coarse problems are proposed. The first two employ a nested coarse space, one obtained by a projection of the fine matrix on the coarse space, and the other obtained by enforcing the equations at a new set of coarse collocation points. The third coarse problem follows the same approach as the first one above, but employs a non-nested coarse space. After presenting our isogeometric OAS preconditioner, we test its performance through a series of numerical tests in two and three dimensions. Since the matrices of collocation schemes are not symmetric, we accelerate our preconditioner with GMRES [32] in order to solve iteratively the resulting collocation problem. We first show the low number of iterations required to solve the OAS preconditioned system (as opposed to the unpreconditioned system), studying how the iteration counts depend on the mesh size and number of subdomains. In particular, from the numerical results one can appreciate the scalability of the two-level OAS method. We also study the method behavior in terms of the polynomial degree p and regularity index k. The iteration counts are very robust with respect to these parameters, even for minimal overlap, a property that is not shared by minimal overlap OAS in Galerkin IGA. We also test the scheme behavior in presence of domain deformation and discontinuous elliptic coefficients, with very satisfactory results. We remark that, although the numerical tests in the current contribution have been implemented on a single-core machine, the proposed OAS preconditioner is also immediately suitable for a parallel implementation on multi-core machines. The present paper is organized as follows. In Section 2 we give a brief presentation of NURBS and isogeometric analysis. In Section 3 we describe the proposed overlapping domain decomposition preconditioner. In Section 4 we present the results of extensive numerical tests in two and three dimensions.
2
2. NURBS-based isogeometric analysis Non-Uniform Rational B-splines are a standard tool for describing and modeling curves and surfaces in computer aided design and computer graphics (see Piegl and Tiller [27] and Rogers [28] for an extensive description of these functions and their properties). In this work, we use NURBS as an analysis tool, as proposed by Hughes et al. [20]. In this section, we present a short description of B-splines, NURBS, the basics of isogeometric analysis and an introduction to the proposed discretization method. In the rest of the paper, we will adopt the following compact notation. Given two real numbers a, b we write a . b, when a ≤ Cb for a generic constant C independent of the knot vectors (defined below), and we write a ≈ b when a . b and b . a. 2.1. B-splines B-splines in the plane are piecewise polynomial curves composed of linear combinations of Bspline basis functions. A knot vector is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve {ξ1 = 0, ..., ξn+p+1 = 1},
(2)
where p is the polynomial degree of the B-spline and n is the number of basis functions (and control points) necessary to describe it. The open interval (ξ1 , ξn+p+1 ) is called a patch. A knot vector is said to be uniform if its knots are uniformly-spaced and non-uniform otherwise. The maximum allowed knot multiplicity is p + 1; a knot vector is said to be open if its first and last knots have multiplicity p + 1. In what follows, we always employ open knot vectors. Basis functions formed from open knot vectors are interpolatory at the ends of the parametric interval Ib := (0, 1) but are not, in general, interpolatory at interior knots. Given a knot vector, univariate B-spline basis functions are defined recursively starting with p = 0 (piecewise constants) ( 1 if ξi ≤ ξ < ξi+1 Ni0 (ξ) = (3) 0 otherwise. For p > 1 : ξ − ξi N p−1 (ξ) + ξi+p+1 − ξ N p−1 (ξ) if ξ ≤ ξ < ξ i i+p+1 p ξi+p+1 − ξi+1 i+1 Ni (ξ) = ξi+p − ξi i 0 otherwise,
(4)
where, in (4), we adopt the convention 0/0 = 0. Thus, the general basis function Nip has support Θi := supp(Nip ) = (ξi , ξi+p+1 ),
i = 1, 2, .., n.
The functions Nip are a partition of unity, see e.g. [14]. In Figure 1 we present an example consisting of n = 9 cubic basis functions generated from the open knot vector ξ = {0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1} If internal knots are not repeated, B-spline basis functions are C p−1 -continuous. If a knot has multiplicity α, the basis is C p−α -continuous at that knot. In particular, when a knot has multiplicity p, the basis is C 0 and interpolates the control point at that location. In the following, we will assume that the maximum knot multiplicity is p so that all considered functions will be (at least) globally continuous. We define the spline space Sbh = span{Nip (ξ), i = 1, . . . , n}. 3
(5)
By means of tensor products, a multi-dimensional B-spline region can be constructed. We discuss here the case of a two-dimensional region, the higher-dimensional case being analogous. b := (0, 1) × (0, 1) be the two-dimensional parametric space. Consider the knot vectors {ξ1 = Let Ω 0, ..., ξn+p+1 = 1} and {η1 = 0, ..., ηm+q+1 = 1}, and an n × m net of control points Ci,j . One dimensional basis functions Nip and Mjq (with i = 1, ..., n and j = 1, ..., m) of degree p and q, b is then defined by respectively, are defined from the knot vectors. The bivariate spline basis on Ω the tensor product construction p,q Bi,j (ξ, η) = Nip (ξ) Mjq (η).
Observe that the two knot vectors {ξ1 = 0, ..., ξn+p+1 = 1} and {η1 = 0, ..., ηm+q+1 = 1} generate a mesh of rectangular elements in the parametric space in a natural way. Analogous to (5), we define p,q Sbh = span{Bi,j (ξ, η), i = 1, . . . , n, j = 1, . . . , m}
(6)
2.2. NURBS maps and spaces
In general, a rational B-spline in Rd is the projection onto d-dimensional physical space of a polynomial B-spline defined in (d+1)-dimensional homogeneous coordinate space. For a complete discussion, see the book by Farin [19] and references therein. In this way, a great variety of geometrical entities can be constructed and, in particular, all conic sections in physical space can be obtained exactly. To obtain a NURBS curve in R2 , we start by introducing the NURBS basis functions of degree p N p (ξ)ωi N p (ξ)ωi Rip (ξ) = Pn i p = i , (7) w(ξ) ˆı=1 Nˆı (ξ)ωˆı P where the denominator w(ξ) = ˆnı=1 Nˆıp (ξ)ωˆı ∈ Sbh is called the weight function. The NURBS curve is then defined by C(ξ) =
n X
Rip (ξ)Ci ,
(8)
i=1
where Ci ∈ R2 are control points. Analogously to B-splines, NURBS basis functions on the two-dimensional parametric space b = (0, 1) × (0, 1) are defined as Ω p,q Ri,j (ξ, η)
p,q p,q Bi,j (ξ, η)ωi,j Bi,j (ξ, η)ωi,j = Pn Pm = , p,q w(ξ, η) ˆı=1 ˆ=1 Bˆı,ˆ (ξ, η)ωˆı,ˆ
(9)
where ωi,j = (Cωi,j )3 and the denominator is the weight function denoted also by w(ξ, η). Observe that the continuity and support of NURBS basis functions are the same as for B-splines. NURBS spaces are the span of the basis functions (9). NURBS regions are defined in terms of the basis functions (9). In particular a single-patch domain Ω is a NURBS region associated with the n × m net of control points Ci,j , and we introduce b → Ω given by the geometrical map F : Ω F(ξ, η) =
n X m X i=1 j=1
4
p,q Ri,j (ξ, η)Ci,j .
(10)
Following the isoparametric approach, the space of NURBS scalar fields on the domain Ω is defined, component by component as the span of the push-forward of the basis functions (9) p,q Nh := span{Ri,j ◦ F−1 , with i = 1, . . . , n; j = 1, . . . , m}.
(11)
As observed for instance in [5], in order to obtain spaces with homogeneous Dirichlet boundary conditions it is sufficient to eliminate the first and last function in each coordinate. In the present paper, in order to keep simple homogeneous Dirichlet boundary conditions but make use of a general structure, we do not exploit such a possibility. Therefore the isogeometric space used in the discretization will be Vb = Sbh and V = Nh living in parameter space and in physical space, respectively and the boundary conditions will be enforced in the discrete problem. The interested reader may find more details on isogeometric analysis as well as many interesting applications in a number of recently published papers mentioned in the introduction. The threedimensional case is analogous and not discussed here. 2.3. Isogeometric collocation In the present section we briefly review the isogeometric collocation method for the present problem (1), introduced in [1], see also [2, 29]. We focus mainly on the bi-dimensional case, the one and three dimensional cases being defined analogously. For simplicity of exposition, we consider directly the case where the collocation points are chosen as the Greville abscissae associated to the knot vectors, see [16]. The choice of collocation points is crucial for the stability and good behavior of the discrete problem; another possible (but less common) choice is given by the Demko abscissae [17]. Let ξ i , i = 1, ..., n, be the Greville abscissae related to the knot vector {ξ1 , ..., ξn+p+1 }: . ξ i = (ξi+1 + ξi+2 + ... + ξi+p )/p .
(12)
Analogously, let η j , j = 1, ..., m, be the Greville abscissae related to the knot vector {η1 , ..., ηm+q+1 }. It is easy to see that ξ 1 = η 1 = 0, ξ n = η m = 1, while all the remaining points are in (0, 1). We define the collocation points τij ∈ Ω by the tensor product structure τij = F(b τij ) ,
τbij = (ξ i , η j ) ∈
b , Ω
for i = 1, ..., n, j = 1, ..., m. Then, the isogeometric collocation problem with Greville abscissae reads: find uh ∈ V such that ( − ∇ · (ρ∇uh )(τij ) = f (τij ) i = 2, ..., n − 1, j = 2, ..., m − 1 , (13) uh (τij ) = 0 (i, j) ∈ ({1, n} × {1, ..., m}) ∪ ({1, ..., n} × {1, m}) . The well-posedness and convergence to the exact solution for problem (13) has been studied, only for the one-dimensional case, in [1]. The theory for the isogeometric collocation method in two and three dimensions is still an open issue; nevertheless, a large array of numerical tests in the literature seem to indicate that the method is stable and convergent in all practical cases. Note that we consider homogeneous Dirichlet boundary conditions only for simplicity; different kind of boundary conditions can be clearly treated in the same fashion.
5
Remark 2.1. In order to ease the presentation of the collocation method we here assume that the knot repetitions are such that the ensuing discrete spline and NURBS spaces are at least C 2 , so that the equations in (13) are well defined. In other words, the polynomial degree p ≥ 3 and the number of repetitions of any knot (excluding the boundary ones) never exceeds p − 2. Should this assumption fail, the method can be easily modified as detailed in [2]; we show an example in Section 4.5. 3. The overlapping Schwarz preconditioner In this section, we construct an isogeometric overlapping Schwarz preconditioner for the isogeometric collocation method (13). 3.1. Subdomains and subspace decomposition We describe first the subdomains and subspace decompositions in 1D and then extend them by tensor products to 2D and 3D. The decomposition is first built for the underlying space of spline functions in parameter space, and then easily extended to the NURBS space in the physical domain. We select from the full set of knots {ξ1 = 0, ..., ξn+p+1 = 1} a subset {ξik , k = 1, . . . , N +1} of (non-repeated) interface knots with ξi1 = 0, ξiN+1 = 1. This subset of interface knots defines a decomposition of the closure of the reference interval Ib = [0, 1] =
[
k=1,..,N
Ibk , with Ibk = (ξik , ξik+1 ),
that we assume to have a similar characteristic diameter H ≈ Hk = diam(Ibk ). The interface knots are thus given by ξik for k = 2, .., N . For each of the interface knots ξik we choose an index 2 ≤ sk ≤ N − 1 (strictly increasing in k) that satisfies sk < ik < sk + p + 1, so that the support of the basis function Nspk intersects both Ibk−1 and Ibk . Note that at least one such sk exists; if it is not unique, any choice can be taken. N
3 i
N
1
0
3 i
1
0
0.2
0.4
ξ
0.6
0.8
0
1
(a) r = 0
0
0.2
0.4
ξ
0.6
0.8
1
(b) r = 1
Figure 1: Cubic basis functions formed from ξ = {0, 0, 0, 0, 1/6, 1/3, 1/2, 2/3, 5/6, 1, 1, 1, 1}. For various r, Vb1 is the span of basis functions drawn with dash-dot and solid lines and Vb2 is the span of basis functions drawn with solid and b In particular, the basis functions in common dashed lines, in two subdomains Ib1 = (0, 1/2) and Ib2 = (1/2, 1) of I. are those drawn with a solid line. Small rectangles at the bottom denote the Greville abscissae related to the knot vector ξ.
6
We then define an overlapping decomposition of Ib in the following way. Let r ∈ N be an integer (called the overlap index) counting the basis functions shared by adjacent subdomains. We define the index sets k = 1, 2, .., N, Θk = {j ∈ N : sk − r ≤ j ≤ sk+1 + r}
with the exception that 1 ≤ j ≤ s2 + r for the index set Θ1 and sN − r ≤ j ≤ n for the index set ΘN and the local spaces Vbk = span{Njp (ξ) : j ∈ Θk }
k = 1, 2, .., N,
(14)
with the analogous exception for Vb1 , VbN . An example is shown in Figure 1. These subspaces form an overlapping decomposition of the spline space Vb ; indeed, 2r + 1 represents the number of basis functions in common (in the univariate case) among “adjacent” local subspaces. For r = 0 we have the minimal overlap consisting of just one common basis function between subspaces. We now define the extended subdomains Ibk′ by Ibk′ =
[
Njp ∈Vbk
supp(Njp ) = (ξsk −r , ξsk+1 +r+p+1 ),
(15)
′ . Let us then introduce a (open) coarse knot vector with the analogous exception for Ib1′ and IbN 0 ξ 0 = {ξ10 = 0, ..., ξN = 1} corresponding to a coarse mesh determined by the subdomains Ibk , c +p+1 i.e. ξ 0 = {ξ1 , ξ2 , . . . , ξp , ξi1 , ξi2 , ξi3 , . . . , ξiN−1 , ξiN , ξiN +1 , ξiN +2 . . . , ξiN +p+1 }, (16)
such that the distance between adjacent distinct knots is of order H, ξ1 = · · · = ξp = ξi1 = 0 and ξiN+1 = ξiN +2 = · · · = ξiN +p+1 = 1 and the associated coarse spline space Vb0 := SbH = span{Ni0,p (ξ), i = 1, ..., Nc }
has the same degree p as Sbh and is thus a subspace of Sbh . Note that, in order to have a smaller coarse space Vb0 , in the definition above we repeat only once all the internal knots (see also Remark 3.1). In more than one dimensions, we just proceed by tensor products. For example, in 2D we define subdomains and overlapping subdomains by Ibk = (ξik , ξik+1 ), b ′ = Ib′ × Ib′ , Ω kl
k
l
Ibl = (ηjl , ηjl+1 ),
b kl = Ibk × Ibl Ω
1 ≤ k ≤ N, 1 ≤ l ≤ M.
N M Moreover, we take the indices {sk }N k=2 associated to the {ξik }k=2 and the analogous indices {sl }l=2 M associated to the {ηjl }l=2 . The local index sets are defined by
Θkl = {(i, j) ∈ N2 : sk − r ≤ i ≤ sk+1 + r, sl − r ≤ j ≤ sl+1 + r} for k = 1, 2, .., N and l = 1, 2, .., M . The local and coarse subspaces now are p,q Vbkl = span{Bi,j (ξ, η) : (i, j) ∈ Θkl }, ◦ p,q
◦ p,q
Vb0 = span{B i,j : B i,j (ξ, η) := Ni0,p (ξ)Mj0,q (η), i = 1, ..., Nc , j = 1, ..., Mc }. 7
The decomposition for the NURBS space V in the physical domain is the trivial extension of the one above. Therefore the local subspaces and the coarse space are, up to the usual modifications for the boundary subdomains, p,q Vkl = span{Ri,j ◦ F−1 : (i, j) ∈ Θkl }, ◦ p,q
◦ p,q
V0 = span{Ri,j ◦ F−1 := (B i,j /w) ◦ F−1 , i = 1, ..., Nc , j = 1, ..., Mc }, where we recall that w is the weight function, see (9). The subdomains in physical space are defined as the image of the subdomains in parameter space with respect to the mapping F b kl ), Ωkl = F(Ω
b ′ ). Ω′kl = F(Ω kl
Remark 3.1. Note that less regular (but still nested) coarse spaces can be considered by including more repetitions of the interior knots. The only requirement for nestedness is that each interior knot is repeated at most as many times as the corresponding knot in the vector defining the fine space. Although we here focus on the more regular case (since it leads to the smaller nested coarse space) the cases with knot repetitions are handled exactly in the same way in what follows, see (18) and (19). 3.2. Matrix form of the preconditioner We are now able to present the preconditioner Schwarz operator, focusing as usual on the bidimensional case. Let the (local and coarse) restriction matrixes Rkl : V → Vkl and R0 : V → V0 be T : V → V , k = 1, .., N , l = 1, .., M defined as the transpose of the natural embedding matrixes Rkl kl T R0 : V0 → V . Let Akl , for all k = 1, .., N , l = 1, .., M , be the square matrix associated to the local collocation problems: find ukl h ∈ Vkl such that ( − ∇ · (ρ∇ukl (i, j) ∈ Θikl , h )(τij ) = f (τij ) ukl h (τij ) = 0
(i, j) ∈ Θ∂kl ,
where the internal and boundary index sets are Θikl = (i, j) ∈ Θkl : (i, j) ∈ / ({1, n} × {1, ..., m}) ∪ ({1, ..., n} × {1, m}) Θ∂kl = Θkl \Θikl .
Note that the above matrix exactly corresponds to T Akl = Rkl A Rkl
(17)
where A is the global collocation matrix associated to the original problem (13). Choice of coarse problem. We propose three choices for the coarse problem. The first two make use of a nested coarse space while the third one makes use of a non-nested coarse space. • First nested choice (A0 ). We simply define, in analogy to (17), the coarse matrix by A0 = R0 A R0T ,
(18)
where R0T is the matrix form of the interpolation operator from the coarse space V0 to the fine space V . Note that the above matrix corresponds to a collocation problem in which each equation in the system is obtained as a weighted sum of the collocation equations of the differential equation on more than one point. 8
• Second nested choice (A′0 ). We consider the collocation matrix ensuing from the coarse basis functions and associated Greville abscissae. We thus define A′0 as the square matrix associated to the coarse collocation problem: find u0h ∈ V0 such that ( − ∇ · (ρ∇u0h )(τij0 ) = f (τij0 ) i = 2, ..., Nc − 1, j = 2, ..., Mc − 1 , u0h (τij0 ) = 0
(i, j) ∈ ({1, Nc } × {1, ..., Mc }) ∪ ({1, ..., Nc } × {1, Mc }), (19) 0 where τij are the Greville abscissae associated to the coarse knot vector (16) by the usual rule (12). In what follows, A′0 denotes the global collocation matrix associated to the coarse problem (19).
• Non-nested choice (A′′0 ). The coarse matrix is defined as in (18), but in this case the coarse space is constituted by the bi-linear or tri-linear piecewise polynomials on the coarse mesh. Since V0 is not a subspace of V , R0T is the matrix form of the L2 projection from V0 to V . Finally, the matrix form of the proposed overlapping Schwarz operator is TOAS = BOAS A, where BOAS is the additive Schwarz preconditioner BOAS = R0T A−1 0 R0 +
N X M X
T Rkl A−1 kl Rkl .
(20)
k=1 l=1
The coarse matrix A0 is to be replaced by the matrix A′0 or A′′0 according to the choice of coarse problem. The three dimensional case is handled analogously. e kl assoMore general multiplicative and hybrid versions as well as projection-like operators T ciated with inexact local solvers could also be used, see [30, 31], but will not be considered here. The use of the BOAS preconditioner (20) for the iterative solution of the discrete system Au = f can also be regarded as replacing it with the preconditioned system TOAS u = g, where g = BOAS f , which can be accelerated by a Krylov subspace method. Remark 3.2. Note that, as is standard in collocation methods, the matrix A is in general not symmetric, even if the original problem is self-adjoint. The same holds for the preconditioned operator; therefore GMRES (instead of CG - conjugate gradient) iterations will need to be used in the solution process. On the other hand, note that the assembly of the local and coarse matrixes Akl , A0 is much cheaper than in Galerkin IGA, exactly as is the case for the global matrix A (see [29]). Moreover such matrixes will be sparser than their Galerkin counterparts. The numerical results presented in the next Section indicate that the convergence rate of our collocation OAS solver is comparable to those of the analogous Galerkin OAS solver studied in our previous work [6], i.e. it is independent of the number of subdomains N (scalability), optimal in the ratio H/h, robust with respect to the IGA polynomial degree p and regularity index k, as well as with respect to the presence of discontinuous elliptic coefficients and domain deformations. Unlike in the IGA Galerkin case, a mathematical analysis of these good convergence rate properties is still missing in the IGA collocation case, as well as a mathematical analysis of the more basic approximation properties of IGA collocation. 9
h 1/8 1/16 1/32 1/64
p=3 29.19 28.56 82.10 327.21
h 1/8 1/16 1/32 1/64
p=3 21.67 84.15 334.15 1334.23
IGA Galerkin p=4 p=5 p=6 240.03 2148.25 20818.97 222.55 1812.58 15132.64 215.00 1700.63 13634.81 381.73 1688.11 13493.15 IGA collocation p=4 p=5 p=6 36.51 54.94 177.74 138.42 205.88 287.38 546.52 810.23 1126.91 2179.25 3228.47 4486.68
p=7 222889.58 129178.94 110155.73 108472.87 p=7 695.15 589.37 1497.23 5955.86
Table 1: Condition numbers computed with respect to the 2-norm of the IGA Galerkin and IGA collocation stiffness matrices on the unit square domain, as a function of the mesh size h and of the spline polynomial degree p. Spline regularity is always maximal, thus k = p − 1.
4. Numerical results In the present section, we test the convergence properties of the isogeometric overlapping Schwarz preconditioner (20) for 2D and 3D model elliptic problems (1) on both parametric (reference square or cube) and physical domains, with zero right-hand side, Dirichlet boundary conditions ex sin(y) in 2D or mixed boundary conditions in 3D. The problem is discretized with isogeometric NURBS spaces with associated mesh size h, polynomial degree p, regularity k, using as a starting basis the Matlab isogeometric library GeoPDEs [15]. The domain is decomposed into N overlapping subdomains of characteristic size H and overlap index r, as described in Section 3.1. The discrete problems are solved by the GMRES method with the isogeometric Schwarz preconditioner (20), with a zero initial guess and as stopping criterion a 10−6 reduction of the relative residual. In the following tests, we study how the convergence rate of the OAS preconditioner depends on h, N, p, k, on the jumps of the elliptic coefficients and on domain deformations. We will mostly consider as collocation points the tensor-product Greville abscissae (12), but in some tests we will also consider the Demko abscissae [17] to show that our OAS preconditioner performs equally well for other choices of collocation points. 4.1. 2D tests: condition number of the collocation IGA stiffness matrix In this test, we study the behavior of the condition number of the IGA collocation stiffness matrix with respect to the mesh size h and the spline polynomial degree p. The domain is the unit square and the regularity of the splines is maximal, thus k = p − 1. The condition number is computed with respect to the 2-norm. The results reported in Table 1 show that the condition number of the IGA collocation stiffness matrix increases as h−2 , as that of the IGA Galerkin stiffness matrix for p = 3. We observe that, for coarse meshes and high p, the condition number of the IGA Galerkin stiffness matrix initially decreases. This unexpected behavior has already been observed in [24]. Up to p = 7, the growth with respect to p of the IGA collocation stiffness matrix seems to be much weaker than that of the IGA Galerkin stiffness matrix, which has an exponential growth, see [24].
10
N 2×2 4×4 8×8 16 × 16 32 × 32
* 12
1/h = 8 A0 A′′0 13 13 -
A′0 13
p = 3, k = 2, r = 0 square domain 1/h = 16 1/h = 32 A0 A′′0 A′0 * A0 A′′0 A′0 * 14 15 15 22 16 19 17 29 16 19 21 29 18 21 25 40 36 17 20 52 53 69 -
* 17 21
1/h = 64 A0 A′′0 A′0 21 25 19 22 25 26 18 22 61 16 19 66 -
* 39 55 73 97 167
1/h = 128 A0 A′′0 A′0 26 34 23 28 34 25 23 25 62 18 21 89 16 19 83
Table 2: 1-level (*) and 2-level OAS preconditioners with nested (A0 ), non-nested (A′′0 ) and second choice (A′0 ) coarse spaces, on the unit square domain: GMRES iteration counts as a function of the number of subdomains N and mesh size 1/h. Fixed p = 3, k = 2, r = 0.
35
21 A’’0, p=4, k=3
20
A’’0, p=3, k=2
19
30 0
GMRES ITERATIONS
GMRES ITERATIONS
A’’ , p=3, k=2 18
A0, p=4, k=3
17 16
A0, p=3, k=2
15 14 13
A’’ , p=4, k=3 0
25
A , p=3, k=2 0
A , p=4, k=3 0
20
15
12 11 10 0
200
400 600 800 NUMBER OF SUBDOMAINS N
1000
10 0
1200
10
20
30
40 50 RATIO H/h
60
70
80
Figure 2: GMRES iteration counts for increasing number of subdomains and fixed H/h = 4 (left panel) and for increasing ratio H/h and fixed N = 2 × 2 (right panel). 2-level Overlapping Schwarz preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse spaces, minimal overlap r = 0, (p = 3, k = 2) and (p = 4, k = 3) spline spaces on the unit square domain.
N 2×2 4×4 8×8 16 × 16 32 × 32
1/h = 8 * A0 A′′0 13 13 12 -
p = 3, k = 2, 1/h = 16 * A0 A′′0 17 14 15 21 16 18 -
r = 0 square domain 1/h = 32 1/h = 64 ′′ * A0 A0 * A0 A′′0 22 16 19 29 20 25 30 18 21 40 22 25 36 17 20 52 19 22 69 16 20 -
1/h = 128 * A0 A′′0 38 27 35 54 28 33 73 24 26 100 18 22 161 16 19
Table 3: 1-level (*) and 2-level OAS preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse spaces, on the unit square domain: GMRES iteration counts as a function of the number of subdomains N and mesh size 1/h. Fixed p = 3, k = 2, r = 0. Demko abscissae are considered as collocation points.
11
N 2×2 4×4 8×8 16 × 16 32 × 32
p = 3, k = 2, r = 0 quarter of ring domain 1/h = 8 1/h = 16 1/h = 32 1/h = 64 * A0 A′′0 * A0 A′′0 * A0 A′′0 * A0 A′′0 16 15 14 21 18 17 28 22 22 37 28 28 28 21 21 39 26 24 55 34 31 53 26 22 76 33 26 104 30 23 -
1/h = 128 * A0 A′′0 50 36 38 73 45 39 109 42 33 171 37 26 245 31 22
Table 4: 1-level (*) and 2-level OAS preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse space in the quarter of ring domain: GMRES iteration counts as a function of the number of subdomains N and mesh size 1/h. Fixed p = 3, k = 2, r = 0.
4.2. 2D tests: OAS scalability in N and optimality in H/h The GMRES iteration counts of the 1-level and 2-level OAS preconditioners on the unit square are reported in Table 2 for fixed p = 3, k = 2, r = 0. Additional results for p = 4, k = 3, r = 0 are plotted in Figure 2. Moreover, Table 3 shows that numerical results obtained by using as collocation points the Demko abscissae [17] are very similar to those in Table 2, indicating that our OAS preconditioners are robust with respect to the choice of collocation points. In what follows, the Greville abscissae are chosen as collocation points. The GMRES iteration counts of the 1-level and 2-level OAS preconditioners on the ring-shaped physical domain for fixed p = 3, k = 2, r = 0 are reported in Table 4. In all the three tables, the iteration counts are reported as function of the number of subdomains N and mesh size 1/h. The results show that the OAS(2) preconditioner is scalable with both nested and non-nested coarse spaces, since moving along the diagonal of the tables the OAS(2) iteration counts are bounded from above by a constant independent of N . This is shown clearly in the left panel of Figure 2. As expected, the coarse problem is essential for scalability, since the results corresponding to the 1-level preconditioner (i.e. without coarse problem) show iteration counts growing with N along the diagonals of the tables. The results also show that the iteration counts depend less than linearly on the ration H/h (right panel of Figure 2 or moving along rows in the tables). These results for varying N and H/h indicate that the proposed collocation OAS preconditioners performs as well as the analogous OAS preconditioners we have previously investigated for Galerkin IGA discretizations (see [6, Sec. 5.1]). 4.3. 2D test: complex spectrum of the OAS preconditioned operator In order to investigate the good convergence rate observed in the previous numerical tests, we report in Figure 3 the spectrum of the 2D OAS preconditioned operator TOAS = BOAS A in the complex plane, considering three cases with increasing number of subdomains N = 16, 64, 144 (first, second, third column, respectively). The increasing ill-conditioning of the unpreconditioned matrix is shown by the spectra of A (top panels), since for increasing N the eigenvalues seem to approach zero and grow unbounded in modulus (for N = 16 and 64 the spectrum of A turns out to be real). On the other hand, the spectrum of the 2-level preconditioner (bottom panels) is contained in a complex region [α, 5] × [−1, 1] with α > 0.5. This fact guarantees the existence of a sequence of complex polynomials with value 1 at the origin and decreasing values on the OAS spectrum, i.e. a fast GMRES convergence, see e.g. [32]. The spectrum of the 1-level preconditioned operator (middle panels) is contained in a similar region but where α seems to approach 0 when N increases, making it impossible to find a sequence of complex polynomials with the same property as before, thus confirming the non-scalability of OAS(1). 12
4
4
2
10
0
2
10
0
0
50 100 150 200 250 min real=9.81e+00 max real=1.39e+03 min im=0.00e+00 max im=0.00e+00
10
300
eig B−1A, OAS(1)
x 10
1 0
eig B
−1 −2
0
0
200 400 600 800 1000 min real=9.88e+00 max real=5.72e+03 min im=0.00e+00 max im=0.00e+00
1200
0
−4
−1
A, OAS(1)
−15
2
0.2 0.1 0 −0.1 −0.2
0.5 1 1.5 2 2.5 3 3.5 4 min real=1.62e−01 max real=4.00e+00 min im=−7.22e−03 max im=7.22e−03
4
x 10
2 0 −2 −4
4.5
0
0.5 1 1.5 2 2.5 3 3.5 4 min real=4.42e−02 max real=4.00e+00 min im=−4.40e−03 max im=4.40e−03
4
eig B
−0.005 −0.01
0
0.5 1 1.5 2 2.5 3 3.5 4 4.5 min real=7.15e−01 max real=4.66e+00 min im=−7.22e−03 max im=7.22e−03
5
x 10
0
4.5
0.5 1 1.5 2 2.5 3 3.5 4 4.5 min real=6.30e−01 max real=4.75e+00 min im=−3.83e−03 max im=3.83e−03
5
0
0 −3
0
−5
0.5 1 1.5 2 2.5 3 3.5 4 min real=1.99e−02 max real=4.00e+00 min im=−3.83e−03 max im=3.83e−03
x 10
−2 −4
4.5
eig B−1A, OAS(2)
eig B−1A, OAS(2)
0
−1
A, OAS(2)
0.005
5
14000
2
−3
0.01
2000 4000 6000 8000 10000 12000 min real=9.90e+00 max real=1.31e+04 min im=−2.05e−01 max im=2.05e−01 −4
eig B−1A, OAS(1)
10
eig A
10
eig A
eig A
10
0.5 1 1.5 2 2.5 3 3.5 4 4.5 min real=6.54e−01 max real=4.73e+00 min im=−4.40e−03 max im=4.40e−03
5
5
x 10
0
−5
0
Figure 3: Spectrum of the collocation unpreconditioned (first row), 1-level (second row) and 2-level (third row) OAS preconditioned operators. 2D quarter-ring domain, N = 16 (first column), 64 (second column), 144 (third column), H/h = 4, p = 3, k = 2.
N k=2 k=3 k=4 k=5 k=6
h = 1/64, N = 4 × 4(H/h = 16), r p=3 p=4 p=5 NP A0 A′′0 NP A0 A′′0 NP A0 215 22 25 307 20 26 414 17 316 23 23 399 22 393 20 -
= 0 square domain p=6 A′′0 NP A0 A′′0 26 467 26 30 28 433 17 26 25 444 21 30 513 20 22 -
NP 655 648 668 473 570
p=7 A0 30 18 18 18 19
A′′0 32 31 26 26 25
Table 5: Unpreconditioned (NP) GMRES and 2-level OAS preconditioners with nested (A0 ) and non-nested (A′′0 ) coarse space, on the unit square domain: GMRES iteration counts as a function of the spline polynomial degree p and the regularity k. Fixed 1/h = 64, N = 4 × 4, H/h = 16, symmetric minimal overlap (r = 0).
13
ρ
NP
OAS(1)
1e-4 1e-2 1e+0 1e+2 1e+4
> 104 1753 496 > 104 > 104
37 48 70 58 62
ρ=1
α
ρ = 10
OAS(2) A0 A′′0 39 39 32 42 30 42 50 42 87 45
Table 6: Central jump test on the unit square domain (left panel). OAS robustness with respect to jump discontinuities in the elliptic coefficient ρ (right table): Iteration counts of the unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with both nested A0 and non-nested A′′0 coarse spaces) preconditioned GMRES. h = 1/64, N = 8 × 8, H/h = 8, overlap r = 0, p = 3, k = 2.
4.4. 2D tests: OAS dependence on p and k The aim of this test is to compare the GMRES iteration counts of the unpreconditioned case and 2-level OAS preconditioners (with nested or non-nested coarse spaces) with respect to the polynomial degree p and spline regularity k. The domain considered is the unit square discretized by a mesh of size h = 1/64, subdivided into N = 4 × 4 subdomains with symmetric minimal overlap r = 0. In this case, if k is odd, we take two common functions instead of one between adjacent univariate subspaces in order to preserve symmetry in terms of basis functions. Table 5 reports the iteration counts, varying p from 3 to 7 and k from 2 to 6. The iteration counts of the unpreconditioned case increase when p increases, while the iteration counts of the 2-level OAS preconditioner seem to have a very robust behavior with respect to k and p. In particular, we note that such behavior is markedly better than its counterpart in [6] for Galerkin IGA discretizations. 4.5. 2D tests: OAS robustness with respect to jumping coefficients In order to investigate the robustness of the OAS preconditioners with respect to jumps discontinuities of the elliptic coefficient ρ, we consider the classical “central jump” test in the unit square decomposed into 8 × 8 subdomains. The elliptic coefficient ρ varies by 8 orders of magnitude (from 10−4 to 104 ) in the 4 × 4 central subdomains, while it is 1 in the surrounding subdomains, see the left figure in Table 6. We fix h = 1/64, H/h = 8, p = 3, k = 2, r = 0. The regularity across the jump interface is reduced to C 0 . Following [2], we enforce the continuity of the weighted normal derivative at the collocation points on the interface layer. Table 6 reports the iteration counts for the unpreconditioned (NP) system, 1-level OAS (OAS(1)) and 2-level OAS (OAS(2)) preconditioned GMRES, with both nested (A0 ) and non-nested (A′′0 ) coarse spaces. The results show clearly the robustness of 1-level OAS and 2-level OAS with non-nested coarse space, since the GMRES iterations remain bounded when ρ jumps. On the contrary, the nested case seems to be more sensitive to high values of the jumps. 4.6. 2D tests: OAS robustness with respect to domain deformations We now study how the OAS preconditioners perform in presence of domain deformations. The four domains in Figure 4 are taken into account, deforming the rectangle [0, 0.5] × [0, 2] (A) into increasingly more curved boomerang-shaped domains (B, C, D). All the domains are discretized by the same mesh of size h = 1/64, subdivided into N = 4× 4 subdomains. The NURBS discretization space is (p = 3, k = 2). Both Greville and Demko collocation points are considered. The results show that both the 1-level and 2-level OAS preconditioners are considerably less sensitive to domain 14
2
2
2
2
1.5
1.5
1.5
1.5
1
1
1
1
0.5
0.5
0.5
0.5
0 0
0.5
(a)
0 0
0.5
1
0 0
(b)
0.5
1
1.5
0 0
0.5
(c)
1
1.5
2
2.5
(d)
Figure 4: Domains for the 2D boomerang test.
domain A B C D
Greville abscissae NP OAS(1) OAS(2) A0 A′′0 194 46 34 32 223 49 39 36 290 50 45 43 492 52 52 51
Demko abscissae NP OAS(1) OAS(2) A0 A′′0 177 45 32 31 214 50 37 36 302 50 45 44 485 52 53 52
Table 7: 2D boomerang test: OAS robustness with respect to domain deformations. Iteration counts of unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with nested A0 and non-nested A′′0 coarse spaces) preconditioned GMRES with Greville and Demko abscissae. h = 1/64, N = 4 × 4, H/h = 16, overlap r = 0, p = 3, k = 2.
deformations than the unpreconditioned solver. Table 7 shows that the GMRES iteration counts increase (going from A to D) by a factor 2.5 for the unpreconditioned solver, while only by a factor 1.5 for 2-level preconditioners. We note that the number of iterations for the 1-level preconditioner seems almost independent of the domain deformation. Nevertheless, since OAS(1) is non-scalable, we expect its iteration counts to increase for increasing subdomains N , while we expect the OAS(2) iteration counts to remain bounded. 4.7. 3D tests: OAS scalability in N The iteration counts of unpreconditioned (NP) GMRES, 1-level OAS (OAS(1)) and 2-level OAS with nested coarse space (OAS(2)) preconditioned GMRES are reported in Table 8 for the reference cubic domain, as a function of the number of subdomains N for fixed subdomain size H/h = 4 (i.e. both h and H are decreasing proportionally as in a scaled speedup test), p = 3, k = 2 and overlap r = 0. These 3D results confirm that 2-level OAS is scalable, since the condition number seems bounded from above by a constant independent of N . 4.8. 3D tests: OAS robustness with respect to domain deformations We finally study how the OAS preconditioners perform in presence of 3D domain deformations. The four domains in Figure 5 are taken into account, deforming the parallelepipedal domain [0, 0.5]× [0, 0.5] × [0, 2] (A) into increasingly more curved boomerang-shaped domains (B, C, D). All the 15
NP 24 33 44 52 64
N 2×2×2 3×3×3 4×4×4 5×5×5 6×6×6
OAS(1) 17 23 29 33 39
OAS(2), A0 17 22 25 26 26
Table 8: OAS scalability in N on the unit cube. Iteration counts of unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with nested A0 coarse space) preconditioned GMRES as a function of the number of subdomains N . Fixed H/h = 4, p = 3, k = 2.
2.5 2 1.5
1.5 1
1
1
0.5
0.5
0.5
0.5
0 2
0 2
0 2
0 2
1.5
1.5
1
1 0.5
0.5
0.5
0.5 0
0
0
(a)
1.5
1
1
(b)
2.5 1.5
1.5
1 0.5
0.5 0 0
0
(c)
2 1.5
1
1
0.5 0
0.5 0
(d)
Figure 5: Domains for the 3D boomerang test.
domains are discretized by the same mesh of size h = 1/12, subdivided into N = 3 × 3 × 3 subdomains. The NURBS discretization space has parameters p = 3, k = 2. Both Greville and Demko collocation points are considered. The results show that, as in the 2D case, both the 1-level and 2-level OAS preconditioners are less sensitive to domain deformations than the unpreconditioned solver. Table 9 shows that the GMRES iteration counts increase by a factor 3 (going from A to D) for the unpreconditioned solver, while only by a factor 1.3 for the 1-level and 2-level preconditioners. References [1] F. Auricchio, L. Beir˜ao da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli. Isogeometric Collocation Methods. Math. Mod. Meth. Appl. Sci., 20 (11), 2075–2107, 2010. [2] F. Auricchio, L. Beir˜ao da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli. Isogeometric collocation for elastostatics and explicit dynamics. Comp. Meth. Appl. Mech. Engrg., 249–252: 2–14, 2012. [3] F. Auricchio, L. Beir˜ao da Veiga, J. Kiendl, C. Lovadina, A. Reali. Locking-free isogeometric collocation methods for spatial Timoshenko rods. Comp. Meth. Appl. Mech. Engrg., 263: 113– 126, 2013.
16
domain A B C D
Greville abscissae NP OAS(1) OAS(2) A0 A′′0 57 32 30 39 68 33 33 33 99 35 37 37 179 38 41 40
Demko abscissae NP OAS(1) OAS(2) A0 A′′0 56 32 30 30 65 33 34 34 91 34 37 36 157 39 41 40
Table 9: 3D boomerang test: OAS robustness with respect to domain deformations. Iteration counts of unpreconditioned (NP) GMRES, 1-level (OAS(1)) and 2-level (OAS(2), with nested A0 and non-nested A′′0 coarse spaces) preconditioned GMRES with Greville and Demko abscissae. h = 1/12, N = 3 × 3 × 3, H/h = 4, overlap r = 0, p = 3, k = 2.
[4] F. Auricchio, F. Calabro, T.J.R. Hughes, A. Reali and G. Sangalli. A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis. Comput. Methods Appl. Mech. Engrg., 249-252: 15–27, 2012. [5] Y. Bazilevs, L. Beir˜ao da Veiga, J.A. Cottrell, T.J.R. Hughes, G. Sangalli. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. Math. Mod. Meth. Appl. Sci., 16, 1–60, 2006. [6] L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. Overlapping Schwarz methods for Isogeometric Analysis. SIAM J. Numer. Anal., 50:1394–1416, 2012. [7] L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. Isogeometric Schwarz preconditioners for linear elasticity systems. Comp. Meth. Appl. Mech. Engrg., 253:439–454, 2013. [8] L. Beir˜ao da Veiga, D. Cho, L.F. Pavarino, S. Scacchi. BDDC preconditioners for Isogeometric Analysis. Math. Mod. Meth. Appl. Sci., 23:1099–1142, 2013. [9] L. Beir˜ao da Veiga, C. Lovadina, A. Reali. Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods. Comp. Meth. Appl. Mech. Engrg., 241–244:38– 51, 2012. [10] L. Beir˜ao da Veiga, L.F. Pavarino, S. Scacchi, O.B. Widlund, S. Zampini. Isogeometric BDDC Preconditioners with Deluxe Scaling. To appear in SIAM J. Sci. Comput. [11] A. Buffa, H. Harbrecht, A. Kunoth, G. Sangalli. BPX-preconditioning for isogeometric analysis. Comp. Meth. Appl. Mech. Engrg., 265:63–70, 2013. [12] N. Collier, D. Pardo, L. Dalcin, M. Paszynsk, V.M. Calo. The cost of continuity: a study of the performance of isogeometric finite elements using direct solvers. Comp. Meth. Appl. Mech. Engrg., 213–216:353–361, 2012. [13] N. Collier, L. Dalcin, D. Pardo and V.M. Calo. The cost of continuity: Performance of iterative solvers on isogeometric finite elements. SIAM J. Sci. Comput., 35, 767–784, 2013. [14] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs. Isogeometric Analysis. Towards integration of CAD and FEA. Wiley, 2009. [15] C. De Falco, A. Reali, R. Vazquez. GeoPDEs: A research tool for Isogeometric Analysis of PDEs. Adv. Engrg. Softw., 42, 1020–1034, 2011. 2010. 17
[16] C. de Boor. A practical guide to splines. Springer, 2001. [17] S. Demko. On the existence of interpolation projectors onto spline spaces, J. of Approx. Theory, 43, 151–156, 1985. [18] M. Dryja, O. B. Widlund. Domain Decomposition Algorithms with Small Overlap, SIAM J. Sci.Comp., 15(3), 604–620, 1994. [19] G.E. Farin. NURBS curves and surfaces: from projective geometry to practical use. A.K. Peters, 1995. [20] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement. Comp. Meth. Appl. Mech. Engrg., 194, 4135– 4195, 2005. [21] T.J.R. Hughes, A. Reali, G. Sangalli. Efficient Quadrature for NURBS-based Isogeometric Analysis. Comp. Meth. Appl. Mech. Engrg., 199, 301–313, 2010. [22] K.P.S. Gahalaut, J.K. Kraus, S.K. Tomar. Multigrid methods for isogeometric discretization. Comp. Meth. Appl. Mech. Engrg., 253:413–425, 2013. [23] K.P.S. Gahalaut, S.K. Tomar, J.K. Kraus. Algebraic multilevel preconditioning in isogeometric analysis: Construction and numerical studies . Comp. Meth. Appl. Mech. Engrg., 266:40–56, 2013. [24] K.P.S. Gahalaut and S.K. Tomar. Condition number estimates for matrices arising in the isogeometric discretizations. Tech. Report RICAM, 23, 2012. [25] H. Gomez, A. Reali, G. Sangalli. Accurate, efficient, and (iso)geometrically flexible collocation methods for phase-field models. Accepted for publication on J. Comp. Phys. [26] S.K. Kleiss, C. Pechstein, B. J¨ uttler, S. Tomar. IETI-Isogeometric Tearing and Interconnecting. Comp. Meth. Appl. Mech. Engrg., 247–248:201–215, 2012. [27] L. Piegl, W. Tiller. The NURBS Book, 2nd Edition. Springer-Verlag, 1997. [28] D.F. Rogers. An Introduction to NURBS With Historical Perspective. Academic Press, 2001. [29] D. Schillinger, J. A. Evans, A Reali, M. A. Scott, T. J.R. Hughes. Isogeometric Collocation: cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations. Comp. Meth. Appl. Mech. Engrg., 267:170–232, 2013. [30] B. F. Smith, P. Bjørstad, W. D. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 1996. [31] A. Toselli, O. B. Widlund. Domain Decomposition Methods: Algorithms and Theory. Computational Mathematics, Vol. 34. Springer-Verlag, Berlin, 2004. [32] L. N. Trefethen, D. Bau. Numerical Linear Algebra. SIAM, 1996.
18
HIGHLIGHTS • We study isogeometric collocation methods which are very recent and promising numerical schemes. • We present an overlapping additive isogeometric collocation discretizations.
Schwarz
method
for
• We show that the overlapping additive Schwarz preconditioner is scalable and very robust.