Numerical and computational aspects of some block-preconditioners for saddle point systems

Highlights

• Element-wise approximation of the Schur complement (EWS) for saddle point matrices.

• Spectral bounds are shown when the Schur complement is symmetric negative definite.

• EWS leads to a numerically and computationally optimal iterative solver.

• Numerical simulations on various computer architectures.

Abstract

Linear systems with two-by-two block matrices are usually preconditioned by block lower- or upper-triangular systems that require an approximation of the related Schur complement. In this work, in the finite element framework, we consider one special such approximation, namely, the element-wise Schur complement. It is sparse and its construction is perfectly parallelizable, making it an appropriate ingredient when building preconditioners for iterative solvers executed on both distributed and shared memory computer architectures. For saddle point matrices with symmetric positive (semi-)definite blocks we show that the Schur complement is spectrally equivalent to the so-constructed approximation and derive spectral equivalence bounds. We also illustrate the quality of the approximation for nonsymmetric problems, where we observe the same good numerical efficiency.

Furthermore, we demonstrate the computational and numerical performance of the corresponding preconditioned iterative solution method on a large scale model benchmark problem originating from the elastic glacial isostatic adjustment model discretized using the finite element method.

Introduction

Performing large scale computer simulations of realistic models is a complex task that consists of several tightly interrelated aspects. Provided that the mathematical model is given, we need to consider the following phases: (1) choose a proper discretization method that guarantees the desired accuracy of the solution, (2) choose numerical solution techniques that are robust and efficient with respect to all involved problem and discretization parameters, (3) choose software implementation that will perform in a fast and scalable manner on current computer architectures. This study addresses the above issues and the interplay between discretization, data structures, numerical solution methods, program implementation and the resulting overall computational efficiency of the simulation.

The problem at hand is the glacial isostatic adjustment (GIA), it is the response of the solid Earth to changes in thickness and extent of ice sheets and glaciers. The GIA process causes subsidence or upliftment of the Earth surface and is currently active both in areas which were previously covered by large ice sheets during the latest glacial period and in currently deglaciating regions of the world. Performing consistent, reliable and fast numerical GIA simulations is becoming an increasingly important facet in the study of how glaciated areas will respond to the present global warming trend.

The Finite element framework has been used for many years in the geophysics community to simulate the GIA process, often using the functionality of general purpose commercial software packages such as ABAQUS, e.g. [1], [2]. However, enhanced numerical GIA simulations are not straightforwardly performed with ABAQUS or other commercial packages since important terms in the continuous GIA model, such as prestress advection, cannot be added directly, leading to the necessity to modify the model in order to be able to use the package. An extensively used implementation was provided by [1], [3] but this was shown to not perform well in the case of material compressibility [4]. Further, in particular for small displacement elasticity formulation, ABAQUS does not handle purely incompressible materials but only nearly incompressible. Also, the currently observed performance of the solvers when utilizing multicore architectures is by far not satisfactory for the target problem, which we illustrate with numerical experiments. Last but not least, the software is available at high cost and usually for limited computer facilities. With the focus on overcoming these problems, the present study is part of an effort aimed at implementing the GIA model for finite elements in its original form, to be described below, and allow to use the full range of values of the involved problem parameters. In addition, we aim at achieving optimal numerical efficiency, expressed in robustness of the solution methods with respect to problem, discretization and method parameters, as well as optimal computational efficiency, expressed in overall work, (nearly) linearly proportional to the number of degrees of freedom. Finally, we aim at good overall scalability when utilizing the nowadays high performance computer platforms. To this end we employ the publicly available software package Deal.II [5] that is flexible enough to provide the desired discretization methods and access to the basic numerical solution methods, but also allows for including user-designed building blocks to be used for facilitating the solution process. As we emphasize large scale and computationally heavy numerical simulations, we consider iterative solution methods to be the viable alternative. The latter automatically entails the need to choose an efficient and robust preconditioning technique, in order to accelerate the convergence of the iterative process.

When discretized, the mathematical model gives rise to linear systems with indefinite matrices of saddle point form, which we solve using state-of-the-art block-structured preconditioners. To achieve the desired overall efficiency, taking into account that the blocks of the system matrix are very large themselves, we use inner iterative solvers. This leads to inner–outer type of solution procedures and raises the necessity to properly choose the inner stopping tolerances so that we balance between computational burden and accurate enough block solvers to ensure the fast convergence of the outer method. In general, such preconditioners require a high quality sparse approximation of a (most often dense) Schur complement matrix.

In this paper we choose one particular sparse approximation of the Schur complement, the so-called element-wise approximation, that is particularly attractive for parallel implementation by construction. We prove the high quality of the approximate Schur complement for the class of symmetric indefinite matrices with a positive semi-definite pivot block. The so-constructed approximation is then used as a part of a block precondioner, applied to the discrete elastic GIA model, for which we test both numerical efficiency (number of iterations) and computational efficiency (time) on multicore CPU systems using the OpenMP (shared memory paradigm) and also the MPI (distributed memory) programming model. We demonstrate that open access numerical linear algebra toolboxes not only enable us to model the target problem in its original formulation but also allow us to implement efficient numerical solution techniques and achieve computational efficiency that outperforms the highly optimized sparse direct solvers in ABAQUS even in serial, and even for two-dimensional problems.

The paper is organized as follows. Section 2 comprises the description of the underlying problem of interest, its discretization and algebraic form. In Section 3 we describe the construction of the approximate Schur complement and analyze its quality when applied to indefinite matrices of saddle point form. We present the experimental setting in Section 4 and the obtained numerical results in Section 5. Some conclusions and an outlook are given in Section 6.