Combining p-multigrid and Multigrid Reduction in Time methods to obtain a scalable solver for Isogeometric Analysis

The use of sequential time integration schemes becomes more and more the bottleneck within large-scale computations due to a stagnation of processor’s clock speeds. In this study, we combine the parallel-in-time Multigrid Reduction in Time method with a p-multigrid method to obtain a scalable solver specifically designed for Isogeometric Analysis. Numerical results obtained for two- and three-dimensional benchmark problems show the overall scalability of the proposed method on modern computer architectures and a significant improvement in terms of CPU timings compared to the use of standard spatial solvers. The use of a p-multigrid method significantly reduces the CPU timings for higher values of the spline degree. The Multigrid Reduction in Time method shows both strong and weak scalability up to 2048 cores. Iteration numbers are independent of the number of time steps, mesh width and spline degree. The use of a p-multigrid method significantly reduces the CPU timings for higher values of the spline degree. The Multigrid Reduction in Time method shows both strong and weak scalability up to 2048 cores. Iteration numbers are independent of the number of time steps, mesh width and spline degree.


Introduction
Since its introduction in [1], Isogeometric Analysis (IgA) has become more and more a viable alternative to the Finite Element Method (FEM). Within IgA, the same building blocks (i.e., B-splines and NURBS) as in Computer Aided Design (CAD) are adopted, which closes the gap between CAD and FEM. In particular, the use of highorder splines results in a highly accurate represention of (curved) geometries and has shown to be advantageous in many applications, like structural mechanics [2], solid and fluid dynamics [3] and shape optimization [4]. Finally, the accuracy per degree of freedom (DOF) compared to FEM is significantly higher with IgA [5].
For time-dependent partial differential equations (PDEs), IgA is typically combined with a traditional time integration scheme within the method of lines. Here, the spatial variables are discretized by adopting IgA, after which the resulting system of ordinary differential equations (ODEs) is integrated in time. However, as with all traditional time integration schemes, the latter part becomes more and more the bottleneck in numerical simulations. When the spatial Vol:.(1234567890) Research Article SN Applied Sciences (2022) 4:163 | https://doi.org/10.1007/s42452-022-05043-7 resolution is increased to improve accuracy, a smaller time step size has to be chosen to ensure stability of the overall method. As clock speeds are no longer increasing, but the core count goes up, the parallelizability of the entire calculation process becomes more and more important to obtain an overall efficient method. As traditional time integration schemes are sequential by nature, new parallel-in-time methods are needed to resolve this problem. The Multigrid Reduction in Time (MGRIT) method [6] is a parallel-in-time algorithm based on multigrid reduction (MGR) techniques [7]. In contrast to space-time methods, in which time is considered as an extra spatial dimension, sequential time stepping is still necessary within MGRIT. Space-time methods have been combined in the literature with IgA [8]. Although very successful, a drawback of such methods is the fact that they are more intrusive on existing codes, while MGRIT just requires a routine to integrate the fully discrete problem from one time instance to the next. Over the years, MGRIT has been studied in detail (see [9]) and applied to a variety of problems, including those arising in optimization [10] and power networks [11].
Recently, the authors applied MGRIT in the context of IgA for the first time in the literature [12]. Here, MGRIT showed convergence for a variety of two-dimensional benchmark problems independent of the mesh width h, the spline degree p of the B-spline basis functions and the number of time steps N t . However, as a standard (diagonally preconditioned) Conjugate Gradient method was adopted for the spatial solves within MGRIT, a significant dependency of the CPU timings on the spline degree was visible. Furthermore, the parallel performance of MGRIT was investigated for a very limited number of cores.
In this paper, we extend the research direction set out in [12] by combining MGRIT with a state-of-the-art p-multigrid method [13] to solve the linear systems arising within MGRIT. CPU timings show that the use of such a solver significantly improves the overall performance of MGRIT, in particular for higher values of p. Furthermore, the parallel performance of the resulting MGRIT method (i.e., strong and weak scalability) is investigated on modern computer architectures, showing significant (and close to ideal) speed-ups up to 2048 cores. This paper is structured as follows: In Sect. 2, a twodimensional model problem and its spatial and temporal discretization are considered. The MGRIT algorithm is then described in Sect. 3. In Sect. 4, the adopted p-multigrid method and its components are presented in more detail. In Sect. 5, numerical results obtained for the considered model problem are analyzed for different values of the mesh width, spline degree and the number of time steps and compared to those obtained in [12]. Furthermore, weak and strong scaling studies of MGRIT when adopting a p-multigrid method are performed. Finally, conclusions are drawn in Sect. 7.

Model problem and discretization
As a model problem, we consider the transient diffusion equation: Here, Ω ⊂ ℝ d denotes a simply connected, Lipschitz domain in d dimensions and f ∈ L 2 (Ω) a source term. The above equation is complemented by initial conditions and homogeneous Dirichlet boundary conditions: First, we discretize Eq. (1) by dividing the time interval [0, T] in N t subintervals of size Δt and applying the -scheme to the temporal derivative, which leads to the following equation to be solved at every time step: for ∈ Ω and k = 0, … , N t . Depending on the choice of , this scheme leads to the backward Euler ( = 1 ), forward Euler ( = 0 ) or second-order accurate Crank-Nicolson ( = 0.5 ) method. By rearranging the terms, the discretized equation can be written as follows: To obtain the variational formulation, let V = H 1 0 (Ω) be the space of functions in the Sobolev space H 1 (Ω) that vanish on the boundary Ω . Equation (5) is multiplied with a test function v ∈ V and the result is then integrated over the domain Ω: where we write u( ) k+1 = u k+1 throughout the remainder of this section to improve readability. Applying integration by parts on the second term on both sides of the equation results in where the boundary integral integral vanishes since v = 0 on Ω . To parameterize the physical domain Ω , a geometry (2) u( , 0) = u 0 ( ), ∈ Ω, Provided that the physical domain Ω is topologically equivalent to the unit square, the geometry can be described by a single geometry function . In case of more complex geometries, a family of functions (m) ( m = 1, … , K ) is defined and we refer to Ω as a multipatch geometry consisting of K patches. For a more detailed description of the spatial discretization in IgA and multipatch constructions, the authors refer to chapter 2 of [1]. At each time step, we express u in Eq. (7) by a linear combination of multivariate B-spline basis functions of order p. Multivariate B-spline basis functions are defined as the tensor product of univariate B-spline basis functions i,p (i = 1, … , N) which are uniquely defined on the parameter domain (0, 1) by an underlying knot vector Ξ = { 1 , 2 , … , N+p , N+p+1 } . Here, N denotes the number of B-spline basis functions and p the spline degree. Based on this knot vector, the basis functions are defined recursively by the Cox-de Boor formula [14], starting from the constant ones: Higher-order B-spline basis functions of order p > 0 are then defined recursively: The resulting B-spline basis functions i,p are non-zero on the interval [ i , i+p+1 ) and possess the partition of unity property. Furthermore, the basis functions are C p−m i -continuous, where m i denotes the multiplicity of knot i . Throughout this paper, we consider a uniform knot otherwise.
vector with knot span size h, where the first and last knot are repeated p + 1 times. As a consequence, the resulting B-spline basis functions are C p−1 continuous and interpolatory at both end points. Figure 1 illustrates both linear (left) and quadratic (right) B-spline basis functions based on such a knot vector. As mentioned previously, the tensor product of univariate B-spline basis functions is adopted for the multidimensional case. Denoting the total number of multivariate B-spline basis functions Φ i,p by N dof , the solution u is thus approximated at each time step as follows: Here, the spline space V h,p is defined, using the inverse of the geometry mapping −1 as pull-back operator, as follows: By setting v = Φ j,p , Eq. (7) can be written as follows: where and denote the mass and stiffness matrix, respectively:

Multigrid Reduction in Time
A traditional (i.e., sequential) time integration scheme would solve Eq. (13) for k = 0, … , N t to obtain the numerical solution at each time instance. In this paper, however, we apply the Multigrid Reduction in Time (MGRIT) method to solve Eq. (13) parallel-in-time. For the ease of where k+1 = ΨΔt . Setting 0 equal to the initial condition u 0 ( ) projected on the spline space V h,p , the time integration method can be written as a linear system of equations: A sequential time integration scheme would correspond to a block-forward solve of this linear system of equations. In this paper, however, we adopt MGRIT to iteratively solve Eq. (17). First, we introduce the two-level MGRIT method, showing similarities with the well-known parareal algorithm [15]. In fact, it can be shown that both methods are equivalent, assuming a specific choice of relaxation [16]. Then, the multilevel variant of MGRIT will be presented in more detail.

Two-level MGRIT method
The two-level MGRIT method combines the use of a cheap coarse-level time integration method with an accurate more expensive fine-level one which can be performed in parallel. That is, the linear system of equations given by Eq. (17) can be solved iteratively by introducing a coarse temporal mesh with time step size Δt C = mΔt F . Here, Δt F coincides with the Δt from the previous sections and m denotes the coarsening factor. Figure 2 illustrates both the fine and coarse temporal discretization. The time instances T 0 , T 1 , … , T N t ∕m are referred to as coarse points (or C-points), while the remaining points are called fine points (or F-points). The description of MGRIT is based on this division of the time instances in both coarse and fine points. By applying a numbering strategy that first numbers the F-points and then the C-points, we can write Eq. (17) as follows: where the matrix can be decomposed as follows: where C and F are identity matrices. The 'ideal' restriction and prolongation operator are then defined as follows: Within MGRIT, the 'ideal' prolongation operator P is typically adopted, while the 'ideal' restriction operator is replaced by R = C . The matrix FF is given by: Note that each solve with Ψ corresponds to a single time step within a coarse interval, which is a completely independent process for each coarse interval and can therefore be performed in parallel. The Schur complement matrix in Eq. (19) is given by: Instead of solving for directly, MGRIT solves for a modified matrix ̃ by replacing the operator (Ψ ) m ≈ Φ , which corresponds to applying a single time step of a coarse time integrator. As a true multigrid method, the building blocks of the MGRIT method consist of relaxation (= fine time stepping) and a coarse grid correction (= coarse time stepping). Relaxation involves solving a linear system of the form  where F and C denote the solution at all F-points and C-points, respectively. Within relaxation, the solution is updated at the F-points based on the given values at the C-points. This time stepping from a coarse point C to all neighbouring fine points is also referred to as F-relaxation [6]. On the other hand, time stepping to a C-point from the previous F-point is referred to as C-relaxation. It should be noted that both types of relaxation are highly parallel and can be combined leading to so-called CF-or FCF-relaxation. Figure 3 illustrates both C-and F-relaxation methods. The coarse grid correction involves solving the linear system of equations which is a sequential procedure by design, but is much cheaper compared to the fine time integration (which can be performed in parallel). Here, the vector C is obtained by applying R on .

Multilevel MGRIT method
The solution procedure described above can be extented to a true multilevel MGRIT method. First, we define a hierarchy of L temporal meshes, where the time step size for the discretization at level l (l = 0, 1, … L) is given by Δt F m l . The total number of levels L is related to the coarsening factor m and the total number of fine steps Δt F by L = log m (N t ) . Let (l) (l) = (l) denote the linear system of equations based on the considered time step size at level l. The MGRIT method can then be written as follows: 1. Apply F-relaxation ( = fine time stepping) on (l) (l) = (l) : 2. Determine the residual at level l and restrict it to level l + 1 using the restriction operator R : 3. Solve Eq. (24) ( = coarse time stepping) to obtain (l+1) : 4. Prolongate the correction using the 'ideal' interpolation operator P and update the solution at level l: Recursive application of this scheme until the coarsest level is reached, leads to a so-called V-cycle. However, as with standard multigrid methods, alternative cycle types (i.e., W-cycles, F-cycles) can be defined. At all levels of the multigrid hierarchy, the operators are obtained by rediscretizing Eq. (1) using a different time step size.

p-multigrid method
Within the MGRIT algorithm, fine time stepping is performed in parallel within each time interval. Assuming a backward Euler time integration scheme, the following linear system of equations is solved within each time interval at every iteration: Throughout this section we will omit the time step index k and write the linear system of equations given by Eq. (25) as follows: In a recent paper by the authors [12], this linear system of equations was solved within MGRIT by means of a (diagonally preconditioned) Conjugate-Gradient method. However, as the condition number of the system matrix increases exponentially in IgA with the spline degree p, the use of standard iterative solvers becomes less efficient for higher values of p. As a consequence, alternative solution techniques have been developed in recent years to overcome this dependency [17].
In this paper, we adopt a p-multigrid method [13] specifically designed for discretizations arising in IgA to solve the linear systems within MGRIT. Within the p-multigrid method, a low-order correction is obtained (at level p = 1 ) to update the solution at the high-order level. Starting from the high-order problem, the following steps are performed [13]:  To approximately solve the residual equation given by Eq.
(29) a single W-cycle of a standard h-multigrid method [18], using canonical prolongation and weighted restriction, is applied. As the level p = 1 corresponds to a loworder Lagrange discretization, an h-multigrid method (using Gauss-Seidel as a smoother) is known to be both efficient and cheap [19]. The resulting p-multigrid adopted throughout this paper is shown in Fig. 4. Note that, we directly restrict the residual at the highorder level to level p = 1 . This aggressive p-coarsening strategy has shown to significantly improve the computational efficiency of the resulting p-multigrid method [20], while maintaining its excellent convergence behavior.
Prolongation and restriction operators based on an L 2 projection are adopted to transfer vectors from the highorder level to the low-order level (and vice versa). These transfer operators have been used extensively in the literature [21][22][23] and are given by: Here, the mass matrix p and transfer matrix p 1 are defined as follows: To prevent the explicit solution of a linear system of equations for each projection step, the consistent mass matrix in both transfer operators is replaced by its lumped counterpart by applying row-sum lumping. Note that, row-sum lumping can be applied within the variational formulation, due to the partition of unity and non-negativity of the B-spline basis functions.
Various choices can be made with respect to the smoother at the high-order level. The use of Gauss-Seidel or (damped) Jacobi as a smoother at level p leads to convergence rates of the resulting multigrid method that depend significantly on the spline degree p [24]. Alternative smoothers have been developed in recent years to overcome this shortcoming [25]. In particular, the use of ILUT factorizations [26] (i.e., as a preconditioner within a preconditioned Richardson iteration) has shown to be very effective in the context of IgA [24] and will therefore be adopted throughout the remainder of this paper. An efficient implementation of ILUT is available in the Eigen library [27]. Once the factorization h,p ≈ h,p h,p is obtained, a single smoothing step is applied as follows: Fig. 4 Illustration of the p-multigrid method [13]. At p = 1 , Gauss-Seidel is adopted as a smoother (filled circle), whereas at the high-order level ILUT is applied (filled triangle). At the coarsest level, a direct solver is applied to solve the residual equation (filled square)

Numerical results
To assess the quality of MGRIT when applied in combination with a p-multigrid method within Isogeometric Analysis, we consider the time-dependent heat equation in two dimensions given by Eq. (1). Figure 5 shows the resulting solution u at different time instances for Ω = [0, 1] 2 . Here, an inhomogeneous Neumann boundary condition is applied at the left boundary. Furthermore, the right-hand side is chosen equal to one and the initial condition is equal to zero. Based on a spatial discretization with B-spline basis functions of order p and a mesh width h, MGRIT is applied to iteratively solve the resulting equation. Both the number of iterations and CPU timings needed to reach convergence will be investigated using both a (diagonally preconditioned) Conjugate Gradient method and the described p-multigrid method. Furthermore, we will investigate the parallel performance of MGRIT on modern computer architectures. The open-source C++ library G+Smo [28] is used to discretize the model problem in space using IgA, while, for the MGRIT algorithm, the parallel-in-time code XBraid, developed at Lawrence Livermore National Lab, is adopted [29]. The MGRIT method is said to have reached convergence if the relative residual (in the L 2 norm) at the end of an iteration is smaller or equal to 10 −10 , unless stated otherwise. As a starting point, we briefly summarize the results obtained in a previous paper of the authors (see [12]). There, numerical results were obtained for the same model problem using different hierarchies (i.e., a V-cycle, F-cycle and two-level method), time integration schemes (i.e., backward Euler, forward Euler and Crank-Nicolson) and domains of interest (see Fig. 6).
In general, it was observed that MGRIT converged in a low number (i.e., 5-10) of iterations, although the number of iterations was slightly higher when V-cycles were adopted instead of F-cycles or a two-level method. Furthermore, the number of iterations was independent of the mesh width h, spline degree of the B-spline basis functions p and the number of time steps N t for all considered hierarchies and domains of interest. As expected from sequential time stepping methods, the use of the implicit backward Euler within MGRIT lead to the most stable time integration method. Finally, CPU timings were obtained for a limited number of processors, showing a strong dependency on the spline degree p when the Conjugate Gradient method was applied as a spatial solver within MGRIT.
In this section, we investigate the effect of using a p-multigrid method for the spatial solves compared to the use of a Conjugate Gradient method. Furthermore, we present numerical results when considering a three dimensional geometry (i.e., the unit cube). Finally, we investigate the weak and strong scaling of MGRIT on modern architectures when applied in the context of IgA. As this research focuses on the spatial solver and scalability, we will restrict ourselves to the backward Euler method and the use of V-cycles within MGRIT.

Iteration numbers
As a first step, we compare the number of MGRIT iterations to reach convergence when a p-multigrid method or a (diagonally preconditioned) Conjugate Gradient method is adopted while keeping all other parameters the same. Table 1 shows the results when the mesh width is kept constant ( h = 2 −6 ) for the unit square and a quarter annulus when adopting V-cycles with a p-multigrid (top) and CG method (bottom), respectively. For both benchmarks and all configurations, the number of iterations needed with MGRIT to reach convergence is independent of the number of time steps N t and spline degree p. Furthermore, the number of MGRIT iterations is identical when adopting a p-multigrid method compared to the use of a Conjugate Gradient method. Table 2 shows the results for different values of the mesh width h when the number of time steps is kept constant ( N t = 100 ) for both benchmarks when adopting V-cycles. The number of MGRIT iterations is independent of the mesh width h and spline degree p. Furthermore, the number of MGRIT iterations is identical when adopting a p-multigrid method compared to the use of a Conjugate Gradient method.    Results when adopting the p-multigrid method have been obtained for a three-dimensional benchmark problem as well. Table 3 shows the number of MGRIT iterations for different values of N t , p and h when the unit cube is considered as geometry. In general, the number of iterations needed to reach convergence is independent of the number of time steps, spline degree and mesh width. Furthermore, the number of iterations are comparable to the ones obtained for the two-dimensional benchmark problems.
Finally, we investigate the influence of the time integration scheme on the number of MGRIT iterations. Table 4 shows the number of MGRIT iterations needed to reach convergence for the forward Euler ( = 0 ) and Crank-Nicolson ( = 0.5 ) method. Results can be compared to the ones obtained with the backward Euler method (see Table 2). For many configurations, MGRIT using forward Euler does not convergence (which is related to the CFL condition), while the Crank-Nicolson method converges for all configurations. A small dependency on h and p is, however, visible. Based on these results, the backward Euler method will be adopted throughout the remainder of this paper. For a more detailed analysis regarding different time integration schemes within MGRIT, the authors refer to [12].
Although the number of MGRIT iterations is identical for all configurations when adopting a p-multigrid or Conjugate Gradient method for solving the linear systems of equations, it is expected that CPU timings will differ significantly. Therefore, focus will lie on CPU timings throughout the remainder of this section.

CPU timings
CPU timings have been obtained when a p-multigrid method or Conjugate Gradient method is adopted for the spatial solves within MGRIT. As in the previous section, we adopt V-cycles, a mesh width of h = 2 −6 and the unit square as our domain of interest. Note that the corresponding iteration numbers can be found in Table 1. The computations are performed on three compute nodes each consisting of an Intel(R) i7-10700 (@ 2.90GHz) Cometlake processor with 8 hardware cores (hyperthreading turned on) and 128GB DDR4 main memory organized in 4 modules of 32GB each. Figure 7 shows the CPU time needed to reach convergence for a varying number of cores, a different number of time steps and different values of p. When the Conjugate Gradient method is adopted for the spatial solves, doubling the number of time steps leads to an increase of the CPU time by a factor of two. Furthermore, it can be observed that the CPU timings significantly increase for higher values of p which is related to the spatial solves required at every time step. As standard iterative solvers (like the Conjugate Gradient method) have a detoriating performance for increasing values of p, more iterations are required to reach convergence for each spatial solve, resulting in higher computational costs of the MGRIT method. When focussing on the number of cores, it can be seen that doubling the number of cores significantly reduces the CPU time needed to reach convergence. More precisely, a reduction of 45-50% can be observed when doubling the number of cores to 6, implying the MGRIT algorithm is highly parallelizable.
As with the use of the Conjugate Gradient method, doubling the number of time steps leads to an increase of the CPU time by a factor of two when a p-multigrid method is adopted. For p = 2 , the use of a p-multigrid method leads to higher CPU timings compared to the use of the Conjugate Gradient method for all values of N t . However, the dependency of the CPU timings on the spline degree is significantly mitigated, which leads to a serious decrease of the CPU timings compared to the use of the Conjugate Gradient method when higher values of p are considered. For example, for N t = 2000 and p = 5 a speed-up of more than a factor of 10 is achieved.
Again, increasing the number of cores from 3 to 6, reduces the CPU time needed to reach convergence by 45-50%. These results show that MGRIT combined with a p-multigrid method leads to an overall more efficient method. Therefore, a larger computer cluster will be considered in the next section to further investigate the scalability of MGRIT (i.e., weak and strong scalability) when combined with a p-multigrid method within IgA.

Scalability
In the previous sections, we applied MGRIT adopting a relatively low number of cores. Here, it was shown that the use of a p-multigrid method significantly reduces the dependency of the CPU timings on the spline degree. In this section, we investigate the scalability of MGRIT (combined with a p-multigrid method) on a modern architecture. More precisely, we will investigate both strong and weak scalability on the Lisa system, one of the nationally used clusters of the Netherlands 1 .

Strong scalability
First, we fix the total problem size and increase the number of cores (i.e., strong scalability). That is, we consider the same benchmark problem as in the previous sections, Here p-multigrid is used for all spatial solves within MGRIT but with a mesh width of h = 2 −6 and a number of time steps N t of 10.000. As before, backward Euler is applied for the time integration and V-cycles are adopted as MGRIT hierarchy. Figure 8 shows the CPU timings needed to reach convergence for a varying number of Intel Xeon Gold 6130 (@ 2.10GHz) processors, where each processor consists of 16 cores. For all values of p, increasing the number of cores leads to significant speed-ups which illustrates the parallizability of the MGRIT method up to 2048 cores. To compare the results with a sequential time integration method, results with a backward Euler method have been added as well. Here, the CPU timings are independent of the number of processors and shown in the most right column for each value of p ('sequential'). Clearly, MGRIT outperforms the sequential algorithm when the number of cores is larger or equal to 128. This behavior has been observed in the literature as well in case of a finite difference discretization for a similar model problem, see [6]. Figure 9 shows the obtained speed-ups as a function of the number of cores for different values of p based on the results presented in Fig. 8. As a comparison, the ideal speed-up has been added, assuming a perfect parallizability of the MGRIT method. Note that, for all values of p, the observed speed-up slightly increases when the number of cores is higher than 256. Furthermore, the obtained speed-ups remain high, even when the number of cores is further increased to 2048, and is independent of the spline degree p.
Strong scalability has been investigated for the threedimensional benchmark problem as well. Figure 10 shows the strong scalability for MGRIT on the unit cube. In general, the obtained results are comparable to the ones obtained in two dimensions, showing significant speedups when increasing the number of cores for all values of p. Results obtained with a sequential time integration method have been added as well, showing comparable results to MGRIT when adopting 512 cores. It should be noted that, compared to the two-dimensional benchmark problem, the CPU timings grow significantly faster for increasing values of p. This is well-known in Isogeometric Analysis [30] and is related to the relatively high number of nonzero entries in three dimensions when considering higher values of p. The use of the (preconditioned) Conjugate Gradient method would even lead to a significant higher growth in CPU timings, as the number of iterations needed to reach convergence for every spatial solve increases excessively in three dimensions when adopting a standard solver. Figure 11 shows the obtained speed-ups for different values of p based on the results presented in Fig. 10. The obtained speed-ups are similar to the ones obtained for the two-dimensional problem but vary slightly more for different values of p. In general, the observed speed-ups remain high, even when the number of cores is increased to 2048.

Weak scalability
As a next step, we consider the unit square as our domain of interest but keep the problem size per processor fixed (i.e. weak scalability). In case of 64 cores, the number of time steps equals 1000 and is adjusted based on the number of cores. Figure 12 shows the CPU time needed to reach convergence for a different number of cores and different values of p. Clearly, the CPU timings remain (more or less) constant when the number of cores is increased, showing the weak scalability of the MGRIT method. Although the CPU timings slightly increase for higher values of p, the strong p-dependency observed with the Conjugate Gradient method is clearly mitigated.

Conclusions
In this paper, we combined MGRIT with a p-multigrid method for discretizations arising in Isogeometric Analysis. Numerical results obtained for a variety of benchmark problems show that the use of a p-multigrid method for all spatial solves within MGRIT results in convergence rates independent of the mesh width h, spline degree p and number of time steps N t . Furthermore, CPU timings depend only mildly on the spline degree p in two dimensions. This is in sharp contrast to standard solvers (e.g. a Conjugate Gradient method), which show a deteriorating performance (in terms of CPU timings) for higher values of p already in two dimensions. Furthermore, the obtained CPU timings when adopting a p-multigrid method are significantly lower for almost all considered configurations. On modern computer architectures, both strong and weak scalability of the resulting MGRIT method have been investigated, showing good scalability up to 2048 cores, illustrating the potential of MGRIT (combined with a p-multigrid method) for time-dependent simulations in IgA.
Within this paper, we restrict ourselves to first-and second-order accurate time integration schemes. As the use of high-order B-spline basis functions significantly reduces the spatial discretization error, the use of alternative (and in particular higher-order) time integration scheme is interesting and will be investigated in future work. Furthermore, we will focus on the application of MGRIT to more challenging benchmark problems, in particular those where IgA has proven to be a viable alternative to FEM.

Author Contributions
All authors contributed to the study conception and design. All authors read and approved the final manuscript.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/