The diploma theses is focused on matrix-free preconditioning of a linear system. It gives a very brief introduction into the area of iterative methods, preconditioning and matrix-free environment. The emphasis is put on a detailed description of a variant of LU factorization which can be computed in a matrix-free manner and on a new technique connected with this factorization for preconditioning by incomplete LU factors in matrix-free environment. Its main features are storage of only one of the two incomplete factors and low memory costs during the computation of the stored factor. The thesis closes with numerical experiments demonstrating the efficiency of the proposed technique.

In this thesis, techniques for efficient utilization of modern computer hardwarefor numerical simulation are considered. In particular, we study techniques for improving the performance of computations using the finite element method. One of the main difficulties in finite-element computations is how to perform the assembly of the system matrix efficiently in parallel, due to its complicated memory access pattern. The challenge lies in the fact that many entries of the matrix are being updated concurrently by several parallel threads. We consider transactional memory, an exotic hardware feature for concurrent update of shared variables, and conduct benchmarks on a prototype multicore processor supporting it. Our experiments show that transactions can both simplify programming and provide good performance for concurrent updates of floating point data. Secondly, we study a matrix-free approach to finite-element computation which avoids the matrix assembly. In addition to removing the need to store the system matrix, matrix-free methods are attractive due to their low memory footprint and therefore better match the architecture of modern processors where memory bandwidth is scarce and compute power is abundant. Motivated by this, we consider matrix-free implementations of high-order finite-element methods for execution on graphics processors, which have seen a revolutionary increase in usage for numerical computations during recent years due to their more efficient architecture. In the implementation, we exploit sum-factorization techniques for efficient evaluation of matrix-vector products, mesh coloring and atomic updates for concurrent updates, and a geometric multigrid algorithm for efficient preconditioning of iterative solvers. Our performance studies show that on the GPU, a matrix-free approach is the method of choice for elements of order two and higher, yielding both a significantly faster execution, and allowing for solution of considerably larger problems. Compared to corresponding CPU implementations executed on comparable multicore processors, the GPU implementation is about twice as fast, suggesting that graphics processors are about twice as power efficient as multicores for computations of this kind.

In this thesis, evaluation of H(div)-conforming finite elements is implemented in a high-performance setting and used to solve the incompressible Navier-Stokes equation, obtaining an exactly point-wise divergence-free velocity field. In particular, the anisotropic Raviart-Thomas tensor-product polynomial space is considered, where the finite element operators are evaluated with quadrature in a matrix-free fashion using sum-factorization on tensor-product meshes. The implementation includes evaluation over elements and faces in two- and three-dimensional space, supporting non-conforming meshes with hanging nodes, and using the contravariant Piola transformation to preserve normal components on element boundaries. In terms of throughput, the implementation achieves up to an order of magnitude faster evaluation of finite element operators compared to a matrix-based evaluation. Correctness is demonstrated with optimal convergence rates for various polynomial degrees, as well as exactly divergence-free solutions for the velocity field.

Distributed parameter systems (DPS) comprise an important class of engineering systems ranging from "traditional" such as tubular reactors, to cutting edge processes such as nano-scale coatings. DPS have been studied extensively and significant advances have been noted, enabling their accurate simulation. To this end a variety of tools have been developed. However, extending these advances for systems design is not a trivial task . Rigorous design and operation policies entail systematic procedures for optimisation and control. These tasks are "upper-level" and utilize existing models and simulators. The higher the accuracy of the underlying models, the more the design procedure benefits. However, employing such models in the context of conventional algorithms may lead to inefficient formulations. The optimisation and control of DPS is a challenging task. These systems are typically discretised over a computational mesh, leading to large-scale problems. Handling the resulting large-scale systems may prove to be an intimidating task and requires special methodologies. Furthermore, it is often the case that the underlying physical phenomena span various temporal and spatial scales, thus complicating the analysis. Stiffness may also potentially be exhibited in the (nonlinear) models of such phenomena. The objective of this work is to design reliable and practical procedures for the optimisation and control of DPS. It has been observed in many systems of engineering interest that although they are described by infinite-dimensional Partial Differential Equations (PDEs) resulting in large discretisation problems, their behaviour has a finite number of significant components , as a result of their dissipative nature. This property has been exploited in various systematic model reduction techniques. Of key importance in this work is the identification of a low-dimensional dominant subspace for the system. This subspace is heuristically found to correspond to part of the eigenspectrum of the system and can therefore be identified efficiently using iterative matrix-free techniques. In this light, only low-dimensional Jacobians and Hessian matrices are involved in the formulation of the proposed algorithms, which are projections of the original matrices onto appropriate low-dimensional subspaces, computed efficiently with directional perturbations.The optimisation algorithm presented employs a 2-step projection scheme, firstly onto the dominant subspace of the system (corresponding to the right-most eigenvalues of the linearised system) and secondly onto the subspace of decision variables. This algorithm is inspired by reduced Hessian Sequential Quadratic Programming methods and therefore locates a local optimum of the nonlinear programming problem given by solving a sequence of reduced quadratic programming (QP) subproblems . This optimisation algorithm is appropriate for systems with a relatively small number of decision variables. Inequality constraints can be accommodated following a penalty-based strategy which aggregates all constraints using an appropriate function , or by employing a partial reduction technique in which only equality constraints are considered for the reduction and the inequalities are linearised and passed on to the QP subproblem . The control algorithm presented is based on the online adaptive construction of low-order linear models used in the context of a linear Model Predictive Control (MPC) algorithm , in which the discrete-time state-space model is recomputed at every sampling time in a receding horizon fashion. Successive linearisation around the current state on the closed-loop trajectory is combined with model reduction, resulting in an efficient procedure for the computation of reduced linearised models, projected onto the dominant subspace of the system. In this case, this subspace corresponds to the eigenvalues of largest magnitude of the discretised dynamical system. Control actions are computed from low-order QP problems solved efficiently online.The optimisation and control algorithms presented may employ input/output simulators (such as commercial packages) extending their use to upper-level tasks. They are also suitable for systems governed by microscopic rules, the equations of which do not exist in closed form. Illustrative case studies are presented, based on tubular reactor models, which exhibit rich parametric behaviour.

A deep scientific understanding of complex physical systems, such as the atmosphere, can be achieved neither by direct measurements nor by numerical simulations alone. Data assimilation is a rigorous procedure to fuse information from a priori knowledge of the system state, the physical laws governing the evolution of the system, and real measurements, all with associated error statistics. Data assimilation produces best (a posteriori) estimates of model states and parameter values, and results in considerably improved computer simulations. The acquisition and use of observations in data assimilation raises several important scientific questions related to optimal sensor network design, quantification of data impact, pruning redundant data, and identifying the most beneficial additional observations. These questions originate in operational data assimilation practice, and have started to attract considerable interest in the recent past. This dissertation advances the state of knowledge in four dimensional variational (4D-Var) - data assimilation by developing, implementing, and validating a novel computational framework for estimating observation impact and for optimizing sensor networks. The framework builds on the powerful methodologies of second-order adjoint modeling and the 4D-Var sensitivity equations. Efficient computational approaches for quantifying the observation impact include matrix free linear algebra algorithms and low-rank approximations of the sensitivities to observations. The sensor network configuration problem is formulated as a meta-optimization problem. Best values for parameters such as sensor location are obtained by optimizing a performance criterion, subject to the constraint posed by the 4D-Var optimization. Tractable computational solutions to this "optimization-constrained" optimization problem are provided. The results of this work can be directly applied to the deployment of intelligent sensors and adaptive observations, as well as to reducing the operating costs of measuring networks, while preserving their ability to capture the essential features of the system under consideration. / Ph. D.