This project is concerned with the development and implementation of a novel preconditioning method for the iterative solution of linear systems that arise in the finite element discretisation of the incompressible Navier-Stokes equations with weakly imposed boundary conditions. In this context we studied an augmented approach where the Schur complement associated with the momentum block of the Navier-Stokes equations has special sparse structure. We follow the standard inf-sup stable method of discretising the Navier-Stokes equations by the Taylor-Hood elements with the Lagrange multiplier constraints discretised using the same order approximation on matching grids. The resulting system of nonlinear equations is solved iteratively by Newton's method. The spectrum of the linearised Oseen's problem, preconditioned by the exact augmentation preconditioner was analysed. Then we developed inexact versions of the preconditioner aimed at achieving optimal scaling of the algorithm in terms of computational resources and wall-clock times. The experimental evaluation of the methodology involve a number of benchmark problems in two and three spatial dimensions. The obtained results demonstrate efficiency, robustness and almost optimal scaling of the solution algorithm with respect to the discrete problem size. We used OOMPH-LIB as a testbed for our experiments. The preconditioning strategies were implemented using OOMPH-LIB's Parallel Block Preconditioning Framework. The initial version of the software was significantly upgraded during the course of this project with newly implemented functionalities to facilitate the rapid development of sophisticated hierarchical design of parallel block preconditioners. Parallel performance analysis of the newly introduced functionalities demonstrate negligible overhead in terms of wall-clock time and the framework demonstrate good weak and strong parallel scaling.

Robust preconditioners on block-triangular and block-factorized form for three types of linear systems of two-by-two block form are studied in this thesis. The first type of linear systems, which are dense, arise from a boundary element type of discretization of crack propagation problems. Numerical experiment show that simple algebraic preconditioning strategies results in iterative schemes that are highly competitive with a direct solution method. The second type of algebraic systems, which are sparse, indefinite and nonsymmetric, arise from a finite element (FE) discretization of the partial differential equations (PDE) that describe (visco)elastic glacial isostatic adjustment (GIA). The Schur complement approximation in the block preconditioners is constructed by assembly of local, exactly computed Schur matrices. The quality of the approximation is verified in numerical experiments. When the block preconditioners for the indefinite problem are combined with an inner iterative scheme preconditioned by a (nearly) optimal multilevel preconditioner, the resulting preconditioner is (nearly) optimal and robust with respect to problem size, material parameters, number of space dimensions, and coefficient jumps. Two approaches to mathematically formulate the PDEs for GIA are compared. In the first approach the equations are formulated in their full complexity, whereas in the second their formulation is confined to the features and restrictions of the employed FE package. Different solution methods for the algebraic problem are used in the two approaches. Analysis and numerical experiments reveal that the first strategy is more accurate and efficient than the latter. The block structure in the third type of algebraic systems is due to a fine-coarse splitting of the unknowns. The inverse of the pivot block is approximated by a sparse matrix which is assembled from local, exactly inverted matrices. Numerical experiments and analysis of the approximation show that it is robust with respect to problem size and coefficient jumps.