ABSTRACT. Isogeometric Analysis [1] has become increasingly popular as an alternative to the Finite Element Method. Solving the resulting linear systems when adopting higher order B-spline basis functions remains a challenging task, as most (standard) iterative methods have a deteriorating preformance for higher values of the approximation order p.

Recently, we succesfully applied p-multigrid methods to discretizations arising in Isogeometric Analysis [2]. In contrast to h-multigrid methods, where each level of the multigrid hierarchy corresponds to a different mesh width h, the p-multigrid hierarchy is constructed based on different approximation orders. The residual equation is then solved at level p=1, enabling the use of efficient solution techniques developed for low-order standard FEM. Numerical results show that the number of iterations needed for convergence is independent of both h and p when the p-multigrid method is enhanced with a smoother based on an Incomplete LU factorization with dual treshold (ILUT). However, a slight dependence on the number of patches has been observed for multipatch geometries.

Since the resulting system matrix has a block structure in case of a multipatch geometry, we consider the use of block ILUT as a smoother. Results indicate that the use of block ILUT can be an efficient alternative to ILUT on multipatch geometries within a heterogeneous HPC framework. Prelimenary results for p-multigrid methods adopting a block ILUT smoother will be presented in this talk.

References [1] T.J.R. Hughes, J.A. Cottrell and Y. Bazilevs, Isogeometric Analysis: CAD, Finite Elements, NURBS, Exact Geometry and Mesh Refinement, Computer Methods in Applied Mechanics and Engineering, 194, 4135 - 4195, 2005

[2] R.Tielen, M. Möller, D. Göddeke and C. Vuik, p-multigrid methods and their comparison to h-multigrid methods within Isogeometric Analysis, Computer Methods in Applied Mechanics and Engineering, 372, 2020

ABSTRACT. Overlapping block smoothers efficiently damp the error contributions from highly oscillatory components within multigrid methods for the Stokes equations but they are computationally expensive. This talk is concentrated on the development and analysis of new block smoothers within geometric multigrid methods for the solution of the Stokes equations discretized on staggered grids. We compare the commonly used Vanka smoother with a non-overlapping variant of the Vanka smoother, the so-called triad-wise smoother. While the latter is computationally cheaper, the convergence depends much more on the implementation than that of the overlapping method. We develop new non-overlapping smoothers, the so-called fourfold triad-wise smoothers, and show their efficiency within multigrid methods to solve the Stokes equations. Local Fourier analysis and a general comparison including the computational cost and the convergence properties of the overlapping and non-overlapping methods will be presented.

ABSTRACT. A multigrid method for solving vector field problems posed in H(curl) is developed and analyzed. The smoothers in the multigrid algorithm are based on nonoverlapping domain decomposition methods and provide the uniform convergence of V-cycle algorithms.

ABSTRACT. We revisit the idea of using polynomial methods to improve a simple iteration for use as a smoother within a multigrid V-cycle for a symmetric positive definite (SPD) system, suggested, e.g., by Adams et al. [1], focusing on the case where the action of one smoother step corresponds to an SPD matrix. A two-level bound going back to Hackbusch is optimized by the lesser-known Chebyshev polynomials of the fourth kind. These polynomials feature heavily in the smoothed aggregation literature, e.g. [2], and actually appear in the article by Adams et al. It does not seem to have been appreciated that these polynomials obey the same three-term recurrence as the Chebyshev polynomials of the first kind (just with a different initial condition), leading to a very simple iteration similar to the usual Chebyshev semi-iteration, also suggested by Adams et al. for use as a smoother. The fourth-kind Chebyshev iteration is actually simpler and avoids the ad-hoc choice of the cut-off of the smallest eigenvalue for the smoother to target. We derive full V-cycle contraction factor bounds for polynomial smoothers depending on the polynomial and an approximation constant involving the single-step smoother, and show that the fourth-kind Chebyshev iteration is quasi-optimal for this bound. While we do not know an elegant closed form for the optimal polynomials for this bound, they can be determined numerically, and implemented with a small modification to the fourth-kind Chebyshev iteration. We show numerical results comparing simple iterations with these two polynomial iterations for a simple model problem.

[1] Mark Adams, Marian Brezina, Jonathan Hu, and Ray Tuminaro. “Parallel multigrid smoothing: polynomial versus Gauss-Seidel”. In: Journal of Computational Physics 188.2 (July 2003), pp. 593–610.

[2] Petr Vaněk and Marian Brezina. “Nearly optimal convergence result for multigrid with aggressive coarsening and polynomial smoothing”. In: Applications of Mathematics 58.4 (Aug. 1, 2013), pp. 369–388.

ABSTRACT. Two-stage Gauss-Seidel preconditioners and smoothers for GMRES-AMG incompressible Navier-Stokes momentum and pressure solvers are described. These are applied to the classical R\"uge-Stuben AMG implemented in Hypre-BoomerAMG and smoothed aggregation SA-AMG in Trilinos-MueLu. Drop tolerances are applied to the sparse lower triangular matrices $L$ arising in the matrix splitting $A = D + L + U$ in order to reduce the number of non-zeros and accelerate the solvers on many-core GPU architectures. The number of small matrix elements are found to increase from fine to coarse levels and thus the efficiency gains are greater for large problems with many levels in the $V$-cycle. Unlike the two-stage algorithm, the solver convergence rate with the standard Gauss-Seidel smoother deteriorates with dropping. The solver compute time is reduced by up to 50\% without a change in the convergence rate. We compare this approach to $L$ and $U$ factors obtained from an ILUTP factorization with dropping applied to reduce the number of non-zeros in $L$. We find that the forward recurrence represented by $L$ is less sensitive to dropping than the backward solve implied by $U$.

ABSTRACT. Greenland and Antarctic ice sheets store most of the fresh water on earth and mass loss from these ice sheets significantly contributes to sea-level rise. The simulation of temperature and velocity of the ice sheets gives rise to large highly nonlinear systems of equations. The solution of the associated tangent problems, arising in Newton’s method, is challenging also because of the strong anisotropy of the meshes. We first consider simulations of the ice velocity of Antarctica and the ice temperature of Greenland. We use one-level Schwarz preconditioners as well as GDSW (Generalized Dryja–Smith–Widlund) type preconditioners from the Trilinos package FROSch (Fast and Robust Schwarz), scaling up to 32 k processor cores (8 k MPI ranks and 4 OpenMP threads) for the finest Antarctica mesh; the corresponding velocity problem contains 566 M degrees of freedom. We then study the coupled velocity and temperature problem for the Greenland ice sheet. To the best of our knowledge, it is the first time that a scalable monolithic two-level preconditioner has been used for this multiphysics problem. We present strong scaling results, up to 4 k MPI ranks, using a monolithic GDSW type preconditioner with decoupled extensions from the FROSch package.

ABSTRACT. In this talk we consider an algebraic monolithic multigrid approach for solving coupled physics systems. In particular, we propose a monolithic method that is generally applicable to volume- coupled problems, that is, problems for which the different physics blocks are defined on the same domain. Examples of volume-coupled problems include thermo-structure-interaction and magnetohydrodynamics (MHD). For such problems, standard aggregation approaches can result in strange or non-standard coarse level operators due to uncorrelated coarsening of different coupled physics blocks. Here we propose that different physics blocks share the same aggregates, ensuring a consistent relationship between the different physics degrees of freedom. Our approach is to use the graph of the sub-block corresponding to a chosen physical system to construct the grid transfer operators for all of the physics sub-blocks. This reuse or cloning of aggregates saves a modest amount of time within the multigrid setup phase and the coupling between different physics blocks is preserved without any extra work. We demonstrate this approach through an implementation in MueLu for various MHD problems and compare the performance with alternative preconditioning methods.

ABSTRACT. Fourth-order differential equations play an important role in many applications in science and engineering, particularly in the modeling of thin plates, beams, and liquid crystals. We present a three-field mixed finite-element formulation for modified forms of the biharmonic problem. This formulation is based on introducing the gradient of the solution as an explicit variable, constrained using a Lagrange multiplier. As a result, the problem is rewritten as a saddle-point system, requiring analysis of the resulting finite-element discretization and construction of optimal linear solvers. Monolithic multigrid solvers for the resulting linear systems are developed. We use standard multigrid V-cycles as preconditioners for FGMRES, with a direct solve at the coarsest level and factor-2 coarsening between all grids. Use of a coupled “star” relaxation scheme is shown to lead to an effective solver in 2D and 3D.

ABSTRACT. Although several numerical approaches for the discretization of constraints in contact mechanics problems are available, mortar methods are widely considered the favorite approach due to their variational consistency, high accuracy, and outstanding robustness [1]. In mortar methods, the contact constraints are usually enforced via a Lagrange multiplier method, such that the solution vector now contains primary and dual unknowns and the linear system to be solved in every Newton step exhibits saddle point structure.

Efficient iterative solvers and preconditioners for saddle point systems are well established in the context of Stokes, Navier-Stokes, or Oseen problems with algebraic multigrid preconditioners taking a leading role. In the context of contact mechanics though, the constraint is imposed locally at the contact interface, posing specific requirements for the aggregation scheme. An exemplary algebraic multigrid method for contact problems can be found in [2].

In this presentation, we will introduce a family of algebraic multigrid preconditioners for saddle point systems arising from mortar contact formulations [3]. As key ingredients, a specialized interface aggregation strategy will be devised to form Lagrange multiplier aggregates along the interface, that are consistent with the aggregates of the underlying volume discretization. Through segregated transfer operators, the saddle point structure is preserved throughout all levels of the AMG hierarchy, such that the contact constraints can be incorporated into the coarse level correction. Numerical experiments will study different variants of the preconditioner and demonstrate the robustness of the proposed methods as well as indicate scalability based on Trilinos/MueLu.

[1] A. Popp, M. Gitterle, M. W. Gee, W. A. Wall: A dual mortar approach for 3D finite deformation contact with consistent linearization, Int. J. Numer. Meth. Engrg., 83(11):1428- 1465, 2010 [2] M. F. Adams: Algebraic multigrid methods for constrained linear systems with applications to contact problems in solid mechanics. Numer. Linear Algebra Appl., 11(2–3):141–153, 2004. [3] T. A. Wiesner, M. Mayr, A. Popp, M. W. Gee, W. A. Wall: Algebraic multigrid methods for saddle point systems arising from mortar contact formulations, submitted for publication.

ABSTRACT. The analysis of the convergence behavior of multilevel methods is typically carried out under the assumption that the problem on the coarsest level is solved exactly. This assumption is, however, not satisfied in practical computations due to the use of an iterative solver on the coarsest level, finite precision arithmetic, or both. In this talk we present an abstract description of multilevel methods which allows inexact solves on the coarsest level and discuss its convergence behavior. In particular, by extending the classic work of Yserentant [H. Yserentant, Acta Numer., 2 (1993), pp. 285-326], we show that even under these weaker assumptions it is still possible to derive a bound on the rate of convergence, which is independent of the number of levels. Further, we consider application of multilevel methods to elliptic partial differential equations and their finite element discretization. We discuss both exact and inexact solvers on the coarsest level. We show that the convergence behavior of the multilevel method with an inexact solver on the coarsest level may depend on the mesh size of the initial triangulation.

ABSTRACT. Multigrid methods are popular for solving linear systems derived from discretizing PDEs. Local Fourier Analysis (LFA) is a technique for investigating and tuning multigrid methods. P-multigrid is popular for high-order or spectral finite element methods, especially on unstructured meshes. In this paper we introduce LFAToolkit.jl, a new Julia package for LFA of high-order finite element methods. LFAToolkit.jl analyzes preconditioning techniques for arbitrary systems of second order PDEs and supports mixed finite element methods. Specifically, in this paper we develop LFA of p-multigrid with arbitrary second-order PDEs using high-order finite element discretizations and investigate tuning Jacobi smoothing for two-grid and multigrid schemes with LFAToolkit.jl. Traditional estimates of the optimal Jacobi smoothing parameter are ill-suited to p-multigrid for high-order elements, but analysis of the two-grid smoothing factor can provide a reasonable underestimate of the optimal Jacobi smoothing parameter for multigrid schemes.

ABSTRACT. Local Fourier analysis is a commonly used tool to assess the quality of geometric multigrid methods for translationally invariant operators. In this this talk we demonstrate how the application of local Fourier analysis can be completely automated and generalized to arbitrary, including non-orthogonal, repetitive structures. To do so we introduce the notion of crystal structures, which allow for a natural representation of almost all translationally invariant operators that are encountered in applications, e.g., discretizations of systems of PDEs, colored domain decomposition approaches and last but not least two- or multigrid hierarchies. We show that based on this definition we are able to fully automate the process of local Fourier analysis both with respect to spatial manipulations of operators as well as the Fourier analysis back-end. The latter is enabled by the introduction and exploitation of a rigorous definition of the underlying crystal structure and a suitable definition of wave functions on this structure. We show that the resulting framework reduces the input required by the user to the specification of the involved operators and their underlying structures in their most natural description and show some applications of the framework.

An open-source Python implementation of aLFA is available at gitlab.com/NilsKintscher/alfa

ABSTRACT. Local Fourier analysis (LFA) serves a prominent role in the prediction of the convergence factor of multigrid methods for discretizations of PDEs. However, few discussions about the implementation of LFA for complicated discretizations of PDEs exist, such as higher-order finite elements, as staggered meshes could lead to complex LFA representations of the grid-transfer operator. In this work, we prove that the LFA representation for d-dimensional PDEs is independent of the placement of degrees of freedom (DoFs). Intuitively speaking, the seeding of the unknowns has no direct effect on the LFA presentation, instead, it serves as a unitary transformation between different representations. Thus, different LFA representations for a given PDE have the same spectrum and norm. Furthermore, we provide a uniform representation in terms of the location of the unknowns, named simple representation, where we allocate all of the DoFs at nodes, resulting in a simple and unified way to compute the symbols of the discrete operators, especially for the grid-transfer operators. This simple representation can contribute to the generalization of the implementation of LFA for different types of discretizations and different problems, especially for higher-order discretization methods.

ABSTRACT. We introduce a randomized algorithm, namely RCHOL, to construct an approximate Cholesky factorization for a given Laplacian matrix (a.k.a., graph Laplacian). The exact Cholesky factorization for the matrix introduces a clique in the associated graph after eliminating every row/column. By randomization, RCHOL keeps a sparse subset of the edges in the clique using a random sampling developed by Spielman and Kyng. We prove RCHOL is breakdown free and apply it to solving large sparse linear systems with symmetric diagonally-dominant matrices. In addition, we parallelize RCHOL based on the nested dissection ordering for shared-memory machines. Numerical experiments demonstrated the robustness and the scalability of RCHOL. For example, our parallel code scaled up to 64 threads on a single node for solving the 3D Poisson equation, discretized with the 7-point stencil on a 1024×1024×1024 grid, or one billion unknowns.

ABSTRACT. To ensure the promised optimal efficiency of algebraic multigrid methods in solving large-scale linear systems of weighted graph Laplacians, the setup phase of the AMG algorithm calls for a reasonable estimate of the error in the current iterative solution. However, this can be expensive to compute. In this talk, we present a nearly-linear scheme to compute one such error estimator by means of a Helmholtz decomposition on the graph. By solving a linear system on the spanning tree and a least-squares problem on the cycle space, a localized error estimator is formed to approximate errors on each edge of the graph. This error estimator can thereafter be used in a variety of $\alpha$AMG algorithms to build adaptive aggregation at a modest cost such that the smooth error is transferred to and represented on the coarse levels approximately. Numerical experiments are presented to demonstrate the effectiveness of path cover adaptive algebraic multigrid (PC-$\alpha$AMG) with the help of a posteriori error estimates for discretized second-order elliptic partial differential equations and real world graph Laplacians problems.

ABSTRACT. We propose a Multilevel Subgraph Preconditioner (MSP) to efficiently solve linear equations in graph Laplacians corresponding to general weighted graphs. The MSP preconditioner combines the ideas of expanded multilevel structure from multigrid (MG) methods and spanning subgraph preconditioners (SSP) from Computational Graph Theory. To start, we expand the original graph based on a multilevel structure to obtain an equivalent expanded graph. Although the expanded graph has a low diameter, which is a favorable property for the SSP, it has negatively weighted edges, which is an unfavorable property for the SSP. We design an algorithm to properly eliminate the negatively weighted edges and theoretically show that the resulting subgraph with positively weighted edges is spectrally equivalent to the expanded graph. Then, we adopt algorithms to find SSP, such as augmented low stretch spanning trees, for the positively weighted expanded graph and, therefore, provide an MSP for solving the original graph Laplacian. MSP is practical to find thanks to the multilevel property and has provable theoretical convergence bounds based on the support theory for preconditioning graphs.