Efficient computation and applications of the Calderón projector

The boundary element method (BEM) is a numerical method for the solution of partial differential equations through the discretisation of associated boundary integral equations.BEM formulations are commonly derived from properties of the Calderón projector, a blocked operator containing four commonly used boundary integral operators. In this thesis, we look in detail at the Calderón projector, derive and analyse a novel use of it to impose a range of boundary conditions, and look at how it can be efficiently computed. Throughout, we present computations made using the open-source software library Bempp, many features of which have been developed as part of this PhD. We derive a method for weakly imposing boundary conditions on BEM, inspired by Nitsche’s method for finite element methods. Formulations for Laplace problems with Dirichlet, Neumann, Robin, and mixed boudary conditions are derived and analysed. For Robin and mixed boundary conditions, the resulting formulations are simpler than standard BEM formulations, and convergence at a similar rate to standard methods is observed. As a more advanced application of this method, we derive a BEM formulation for Laplace’s equation with Signorini contact conditions. Using the weak imposition framework allows us to naturally impose this more complex boundary condition; the ability to do this is a significant advantage of this work. These formulations are derived and analysed, and numerical results are presented. Using properties of the Calderón projector, methods of operator preconditioning for BEM can be derived. These formulations involve the product of boundary operators. We present the details of a discrete operator algebra that allows the easy calculation of these products on the discrete level. This operator algebra allows for the easy implementation of various formulations of Helmholtz and Maxwell problems, including regularised combined field formulations that are immune to ill-conditioning near eigenvalues that are an issue for other formulations. We conclude this thesis by looking at weakly imposing Dirichlet and mixed Dirichlet–Neumann boundary condition on the Helmholtz equation. The theory for Laplace problems is extended to apply to Helmholtz problems, and an application to wave scattering from multiple scatterers is presented

3D simulation of magneto-mechanical coupling in MRI scanners using high order FEM and POD

Magnetic Resonance Imaging (MRI) scanners have become an essential tool in the medi-cal industry due to their ability to produce high resolution images of the human body. To generate an image of the body, MRI scanners combine strong static magnetic fields with transient gradient magnetic fields. The interaction of these magnetic fields with the con-ducting components present in superconducting MRI scanners gives rise to an important problem in the design of new MRI scanners. The transient magnetic fields give rise to the appearance of eddy currents in conducting components. These eddy currents, in turn, result in electromagnetic stresses, which cause the conducting components to deform and vibrate. The vibrations are undesirable as they lead to a deterioration in image quality (with image artefacts) and to the generation of noise, which can cause patient discomfort. The eddy currents, in addition, lead to heat being dissipated and deposited into the cryo-stat, which is filled with helium in order to maintain the coils in a superconducting state. This deposition of heat can cause helium boil off and potentially result in a costly magnet quench. Understanding the mechanisms involved in the generation of these vibrations and the heat being deposited into the cryostat are, therefore, key for a successful MRI scanner design. This involves the solution of a coupled magneto-mechanical problem, which is the focus of this work.In this thesis, a new computational methodology for the solution of three-dimensional (3D) magneto-mechanical coupled problems with application to MRI scanner design is presented. To achieve this, first an accurate mathematical description of the magneto-mechanical coupling is presented, which is based on a Lagrangian formulation and the assumption of small displacements. Then, the problem is linearised using an AC-DC splitting of the fields, and a variational formulation for the solution of the linearised prob-lem in a time-harmonic setting is presented. The problem is then discretised using high order finite elements, where a combination of hierarchical H1 and H(curl) basis func-tions is used. An efficient staggered algorithm for the solution of the coupled system is proposed, which combines the DC and AC stages and makes use of preconditioned iter-ative solvers when appropriate. This finite element methodology is then applied to a set of challenging academic and industrially relevant problems in order to demonstrate its accuracy and efficiency.This finite element methodology results in the accurate and efficient solution of the magneto-mechanical problem of interest. However, in the design stage of a new MRI scanner, this coupled problem must be solved repeatedly for varying model parameters such as frequency or material properties. Thus, even if an efficient finite element solver is available for the solution of the coupled problem, the need for these repeated simulations result in a bottleneck in terms of computational cost, which leads to an increase in design time and its associated financial implications. Therefore, in order to optimise this process, the application of Reduced Order Modelling (ROM) techniques is considered. A ROM based on the Proper Orthogonal Decomposition (POD) method is presented and applied to a series of challenging MRI configurations. The accuracy and efficiency of this ROM is demonstrated by performing comparisons against the full order or high fidelity finite element software, showing great performance in terms of computational speed-up, which has major benefits in the optimisation of the design process of new MRI scanners

Accelerated Calderón Preconditioning for Electromagnetic Scattering by Multiple Absorbing Dielectric Objects

We consider electromagnetic scattering by multiple absorbing dielectric objects using the PMCHWT boundary integral equation formulation. Galerkin discretisation of this formulation leads to ill-conditioned linear systems, and Calderón preconditioning, an operator-based approach, can be used to remedy this. To obtain a stable discretisation of the operator products that arise in this approach, the use of a dual mesh defined on a barycentrically refined grid needs to be considered, increasing memory consumption. Furthermore, to capture the oscillatory solution of the electromagnetic waves, the mesh needs to be refined with respect to frequency, making the simulation of high-frequency problems very expensive. This thesis presents two complementary approaches to minimising memory cost and computation time (for assembly and solution): modification of the preconditioning operator, and a bi-parametric implementation. The former aims to minimise the number of operators used in the preconditioner to reduce the additional matrix-vector products performed, and the memory cost, while still maintaining a sufficient preconditioning effect. The latter uses two distinct sets of parameters during assembly, to minimise assembly and solution time as well as memory. The operator is assembled with a more expensive set of parameters to obtain an accurate solution. The preconditioner, which is discretised using the expensive dual basis functions, is assembled with a cheaper set of parameters. The two approaches are explained in the context of a series of model problems, then applied to realistic ice crystal configurations found in cirrus clouds. They are shown to deliver a reduction of 99% in memory cost and at least 80% in computation time, for the highest frequency considered. The accelerated formulations have been used at the Met Office to create a new database of the scattering properties of atmospheric ice crystals for future numerical weather prediction. A brief description of that work is also presented in the thesis

Preconditioning for hyperelasticity-based mesh optimisation

A robust mesh optimisation method is presented that directly enforces the resulting deformation to be orientation preserving. Motivated by aspects from mathematical elasticity, the energy functional of the mesh deformation can be related to a stored energy functional of a hyperelastic material. Formulating the functional in the principal invariants of the deformation gradient allows fine grained control over the resulting deformation. Solution techniques for the arising nonconvex and highly nonlinear system are presented. As existing preconditioners are not sufficient, a PDE-based preconditioner is developed

Template-Based Image Reconstruction from Sparse Tomographic Data

Funder: University of CambridgeAbstract: We propose a variational regularisation approach for the problem of template-based image reconstruction from indirect, noisy measurements as given, for instance, in X-ray computed tomography. An image is reconstructed from such measurements by deforming a given template image. The image registration is directly incorporated into the variational regularisation approach in the form of a partial differential equation that models the registration as either mass- or intensity-preserving transport from the template to the unknown reconstruction. We provide theoretical results for the proposed variational regularisation for both cases. In particular, we prove existence of a minimiser, stability with respect to the data, and convergence for vanishing noise when either of the abovementioned equations is imposed and more general distance functions are used. Numerically, we solve the problem by extending existing Lagrangian methods and propose a multilevel approach that is applicable whenever a suitable downsampling procedure for the operator and the measured data can be provided. Finally, we demonstrate the performance of our method for template-based image reconstruction from highly undersampled and noisy Radon transform data. We compare results for mass- and intensity-preserving image registration, various regularisation functionals, and different distance functions. Our results show that very reasonable reconstructions can be obtained when only few measurements are available and demonstrate that the use of a normalised cross correlation-based distance is advantageous when the image intensities between the template and the unknown image differ substantially

Discontinuous Galerkin Methods for the Linear Boltzmann Transport Equation

Radiation transport is an area of applied physics that is concerned with the propagation and distribution of radiative particle species such as photons and electrons within a material medium. Deterministic models of radiation transport are used in a wide range of problems including radiotherapy treatment planning, nuclear reactor design and astrophysics. The central object in many such models is the (linear) Boltzmann transport equation, a high-dimensional partial integro-differential equation describing the absorption, scattering and emission of radiation. In this thesis, we present high-order discontinuous Galerkin finite element discretisations of the time-independent linear Boltzmann transport equation in the spatial, angular and energetic domains. Efficient implementations of the angular and energetic components of the scheme are derived, and the resulting method is shown to converge with optimal convergence rates through a number of numerical examples. The assembly of the spatial scheme on general polytopic meshes is discussed in more detail, and an assembly algorithm based on employing quadrature-free integration is introduced. The quadrature-free assembly algorithm is benchmarked against a standard quadrature-based approach, and an analysis of the algorithm applied to a more general class of discontinuous Galerkin discretisations is performed. In view of developing efficient linear solvers for the system of equations resulting from our discontinuous Galerkin discretisation, we exploit the variational structure of the scheme to prove convergence results and derive a posteriori solver error estimates for a family of iterative solvers. These a posteriori solver error estimators can be used alongside standard implementations of the generalised minimal residual method to guarantee that the linear solver error between the exact and approximate finite element solutions (measured in a problem-specific norm) is below a user-specified tolerance. We discuss a family of transport-based preconditioners, and our linear solver convergence results are benchmarked through a family of numerical examples

Tensor approximation methods for stochastic problems

Spektrale stochastische Methoden haben sich als effizientes Werkzeug zur Modellierung von Systemen mit Unsicherheiten etabliert. Der Vorteil dieser Methoden ist, dass sie nicht nur Statistiken liefern, sondern auch eine direkte Darstellung der Lösung als sogenanntes Surrogatmodell. Besonders attraktiv für elliptische stochastische partielle Differentialgleichungen (SPDGln) ist das stochastische Galerkin Verfahren, da in diesem wesentliche Eigenschaften des Differentialoperators erhalten bleiben. Ein Nachteil der Methode ist jedoch, dass enorme Mengen an Speicherplatz benötigt werden, da die Lösung in einem Tensorprodukt der räumlichen und stochastischen Ansatzräume liegt. Bisher wurden verschiedene Ansätze erprobt, um diese Anforderung zu verringern. Hierzu zählen Modellreduktionstechniken, Unterraumiterationen, um den Lösungsraum auf einen beherrschbaren Unterraum einzuschränken, oder Methoden, welche die Lösung schrittweise aus Rang-1 Produkten aufzubauen. In der vorliegenden Arbeit werden Bestapproximationen der Lösungen linearer SPDGln als Niedrig-Rang-Darstellungen gesucht. Dies wird dadurch erreicht, dass Tensordarstellungen sowohl für die Eingangsdaten als auch für die Lösung verwendet und während des ganzen iterativen Lösungsprozesses beibehalten werden. Da diese Darstellungen weitere Näherungen während des Lösungsprozesses erfordern, ist es wesentlich die Konvergenz der Lösung genau zu überwachen. Ferner müssen Besonderheiten der Präkonditionierung der diskreten Systeme und der Stagnation der iterativen Verfahren beachtet werden. Mit dem Ziel der praktischen Anwendbarkeit als einem wesentlichen Bestandteil dieser Arbeit wurde großer Wert auf eine detaillierte Beschreibung der Implementierungstechniken gelegt.Spectral stochastic methods have gained wide acceptance as a tool for efficient modelling of uncertain stochastic systems. The advantage of those methods is that they provide not only statistics, but give a direct representation of the measure of the solution as a so-called surrogate model, which can be used for very fast sampling. Especially attractive for elliptic stochastic partial differential equations (SPDEs) is the stochastic Galerkin method, since it preserves essential properties of the differential operator. One drawback of the method is, however, that it requires huge amounts of memory, as the solution is represented in a tensor product space of spatial and stochastic basis functions. Different approaches have been investigated to reduce the memory requirements, for example, model reduction techniques using subspace iterations to reduce the approximation space or methods of approximating the solution from successive rank-1 updates. In the present thesis best approximations to the solutions of linear elliptic SPDEs are constructed in low-rank tensor representations. By using tensor formats for all random quantities, the best subsets for representing the solution are computed “on the fly” during the entire process of solving the SPDE. As those representations require additional approximations during the solution process it is essential to control the convergence of the solution. Furthermore, special issues with preconditioning of the discrete system and stagnation of the iterative methods need adequate treatment. Since one goal of this work was practical usability, special emphasis has been given to implementation techniques and their description in the necessary detail

Tensor approximation methods for stochastic problems

Spektrale stochastische Methoden haben sich als effizientes Werkzeug zur Modellierung von Systemen mit Unsicherheiten etabliert. Der Vorteil dieser Methoden ist, dass sie nicht nur Statistiken liefern, sondern auch eine direkte Darstellung der Lösung als sogenanntes Surrogatmodell. Besonders attraktiv für elliptische stochastische partielle Differentialgleichungen (SPDGln) ist das stochastische Galerkin Verfahren, da in diesem wesentliche Eigenschaften des Differentialoperators erhalten bleiben. Ein Nachteil der Methode ist jedoch, dass enorme Mengen an Speicherplatz benötigt werden, da die Lösung in einem Tensorprodukt der räumlichen und stochastischen Ansatzräume liegt. Bisher wurden verschiedene Ansätze erprobt, um diese Anforderung zu verringern. Hierzu zählen Modellreduktionstechniken, Unterraumiterationen, um den Lösungsraum auf einen beherrschbaren Unterraum einzuschränken, oder Methoden, welche die Lösung schrittweise aus Rang-1 Produkten aufzubauen. In der vorliegenden Arbeit werden Bestapproximationen der Lösungen linearer SPDGln als Niedrig-Rang-Darstellungen gesucht. Dies wird dadurch erreicht, dass Tensordarstellungen sowohl für die Eingangsdaten als auch für die Lösung verwendet und während des ganzen iterativen Lösungsprozesses beibehalten werden. Da diese Darstellungen weitere Näherungen während des Lösungsprozesses erfordern, ist es wesentlich die Konvergenz der Lösung genau zu überwachen. Ferner müssen Besonderheiten der Präkonditionierung der diskreten Systeme und der Stagnation der iterativen Verfahren beachtet werden. Mit dem Ziel der praktischen Anwendbarkeit als einem wesentlichen Bestandteil dieser Arbeit wurde großer Wert auf eine detaillierte Beschreibung der Implementierungstechniken gelegt.Spectral stochastic methods have gained wide acceptance as a tool for efficient modelling of uncertain stochastic systems. The advantage of those methods is that they provide not only statistics, but give a direct representation of the measure of the solution as a so-called surrogate model, which can be used for very fast sampling. Especially attractive for elliptic stochastic partial differential equations (SPDEs) is the stochastic Galerkin method, since it preserves essential properties of the differential operator. One drawback of the method is, however, that it requires huge amounts of memory, as the solution is represented in a tensor product space of spatial and stochastic basis functions. Different approaches have been investigated to reduce the memory requirements, for example, model reduction techniques using subspace iterations to reduce the approximation space or methods of approximating the solution from successive rank-1 updates. In the present thesis best approximations to the solutions of linear elliptic SPDEs are constructed in low-rank tensor representations. By using tensor formats for all random quantities, the best subsets for representing the solution are computed “on the fly” during the entire process of solving the SPDE. As those representations require additional approximations during the solution process it is essential to control the convergence of the solution. Furthermore, special issues with preconditioning of the discrete system and stagnation of the iterative methods need adequate treatment. Since one goal of this work was practical usability, special emphasis has been given to implementation techniques and their description in the necessary detail

Discontinuous Galerkin Methods for the Linear Boltzmann Transport Equation

Radiation transport is an area of applied physics that is concerned with the propagation and distribution of radiative particle species such as photons and electrons within a material medium. Deterministic models of radiation transport are used in a wide range of problems including radiotherapy treatment planning, nuclear reactor design and astrophysics. The central object in many such models is the (linear) Boltzmann transport equation, a high-dimensional partial integro-differential equation describing the absorption, scattering and emission of radiation. In this thesis, we present high-order discontinuous Galerkin finite element discretisations of the time-independent linear Boltzmann transport equation in the spatial, angular and energetic domains. Efficient implementations of the angular and energetic components of the scheme are derived, and the resulting method is shown to converge with optimal convergence rates through a number of numerical examples. The assembly of the spatial scheme on general polytopic meshes is discussed in more detail, and an assembly algorithm based on employing quadrature-free integration is introduced. The quadrature-free assembly algorithm is benchmarked against a standard quadrature-based approach, and an analysis of the algorithm applied to a more general class of discontinuous Galerkin discretisations is performed. In view of developing efficient linear solvers for the system of equations resulting from our discontinuous Galerkin discretisation, we exploit the variational structure of the scheme to prove convergence results and derive a posteriori solver error estimates for a family of iterative solvers. These a posteriori solver error estimators can be used alongside standard implementations of the generalised minimal residual method to guarantee that the linear solver error between the exact and approximate finite element solutions (measured in a problem-specific norm) is below a user-specified tolerance. We discuss a family of transport-based preconditioners, and our linear solver convergence results are benchmarked through a family of numerical examples

Compatible finite element methods for geophysical fluid dynamics

This article surveys research on the application of compatible finite element methods to large scale atmosphere and ocean simulation. Compatible finite element methods extend Arakawa's C-grid finite difference scheme to the finite element world. They are constructed from a discrete de Rham complex, which is a sequence of finite element spaces which are linked by the operators of differential calculus. The use of discrete de Rham complexes to solve partial differential equations is well established, but in this article we focus on the specifics of dynamical cores for simulating weather, oceans and climate. The most important consequence of the discrete de Rham complex is the Hodge-Helmholtz decomposition, which has been used to exclude the possibility of several types of spurious oscillations from linear equations of geophysical flow. This means that compatible finite element spaces provide a useful framework for building dynamical cores. In this article we introduce the main concepts of compatible finite element spaces, and discuss their wave propagation properties. We survey some methods for discretising the transport terms that arise in dynamical core equation systems, and provide some example discretisations, briefly discussing their iterative solution. Then we focus on the recent use of compatible finite element spaces in designing structure preserving methods, surveying variational discretisations, Poisson bracket discretisations, and consistent vorticity transport.Comment: correction of some typo