Research Classification
Research Interests
Relevant Thesis-Based Degree Programs
Affiliations to Research Centres, Institutes & Clusters
Graduate Student Supervision
Doctoral Student Supervision
Dissertations completed in 2010 or later are listed below. Please note that there is a 6-12 month delay to add the latest dissertations.
Analysis and preconditioning of double saddle-point systems (2022)
This thesis deals with the mathematical analysis and numerical solution of double saddle-point systems.We derive bounds on the eigenvalues of a generic form of double saddle-point matrices with a positive definite leading block. The bounds are expressed in terms of extremal eigenvalues and singular values of the associated block matrices. Inertia and algebraic multiplicity of eigenvalues are considered as well. The analysis includes bounds for preconditioned matrices based on block diagonal preconditioners using Schur complements, and it is shown that in this case the eigenvalues are clustered within a few intervals bounded away from zero. Analysis for approximations of Schur complements is included. Some numerical observations validate our analytical findings.We also derive bounds on the eigenvalues of (classical) saddle-point matrices with singular leading blocks. The technique of proof is based on augmentation. Our bounds depend on the principal angles between the ranges or kernels of the matrix blocks. We use these analyses to derive a preconditioner for saddle-point systems with singular leading blocks. Our preconditioning approach is based on augmenting the leading block and using Schur complements of the augmented system. We show that the resulting preconditioned operator has four distinct eigenvalues, and numerical experiments validate the effectiveness of our approach.We then extend the preconditioner for saddle-point systems with a singular leading block to deal with double saddle-point systems with a singular leading block. The preconditioner is based on augmenting the leading block by a null matrix of one of the off-diagonal blocks, and using the Schur complements of the augmented system.
View record
Optimal mapping with topology change for animation and geometry problems (2019)
This thesis investigates several aspects of the optimal mapping problem, finding a bijective function between two shapes which minimizes some metric, and its many applications in computer graphics such as finding planar embeddings of curved surface meshes, image warping, mesh deformation, elastoplastic simulation, and more. We highlight in particular problems where differences in topology between the shapes—which necessitate discontinuities in the mapping—play a crucial role, e.g. due to animation decisions or physical fracture. Our first project presents an approach to extend discrete variational shape interpolation to create keyframe animation involving extreme deformation, topology change and dynamics. To construct correspondence between keyframe shapes as well as satisfying feature matching constraints, we introduce a consistent multimesh structure that is able to resolve topological in-equivalence between shapes and a Comesh Optimization algorithm that optimizes our multimesh for both intra and inter-mesh quality. To further speed up the total solving time, we then develop three new improvements to the state of the art nonlinear optimization techniques that are frequently used in this field, including a barrier-aware line-search filter that cures blocked descent steps due to element barrier terms; an energy proxy model that adaptively blends the Sobolev gradient and L-BFGS descent to gain the advantages of both; and a characteristic gradient norm providing a robust and largely mesh- and energy independent convergence criterion. Finally, we demonstrate another related interesting graphics application, brittle fracture simulation. In particular, we investigate whether we can sidestep the volumetric optimal mapping problem for solving elastic deformation and instead use a boundary element discretization which only requires a surface mesh. We will show how the computational cost of such problems can be alleviated using a boundary integral formulation and kernel-independent Fast Multipole Method. By combining with an explicit surface tracking framework, we further avoid the expensive volumetric mesh construction and maintenance during fracture propagation and shape topology change.
View record
Preconditioners for incompressible magnetohydrodynamics (2019)
The main goal of this thesis is to design efficient numerical solutions to incompressible magnetohydrodynamics (MHD) problems, with focus on the solution of the large and sparse linear systems that arise. The MHD model couples the Navier-Stokes equations that govern fluid dynamics and Maxwell's equations which govern the electromagnetic effects.We consider a mixed finite element discretization of an MHD model problem. Upon discretization and linearization, a large block 4-by-4 nonsymmetric linear system needs to be (repeatedly) solved. One of the principal challenges is the presence of a skew-symmetric term that couples the fluid velocity with the magnetic field.We propose two distinct preconditioning techniques. The first approach relies on utilizing and combining effective solvers for the mixed Maxwell and the Navier-Stokes sub-problems. The second approach is based on algebraic approximations of the inverse of the matrix of the linear system. Both approaches exploit the block structure of the discretized MHD problem.We perform a spectral analysis for ideal versions of the proposed preconditioners, and develop and test practical versions.Large-scale numerical results for linear systems of dimensions up to 10⁷ in two and three dimensions validate the effectiveness of our techniques.We also explore the use of the Conjugate Gradient (CG) method for saddle-point problems with an algebraic structure similar to the time-harmonic Maxwell problem. Specifically, we show that for a nonsingular saddle-point matrix with a maximally rank-deficient leading block, there are two sufficient conditions that allow for CG to be used.An important part of the contributions of this thesis is the development of numerical software that utilizes state-of-the-art software packages. This software is highly modular, robust, and flexible.
View record
Tackling small-scale fluid features in large domain for computer graphics (2017)
Turbulent gaseous phenomena, often visually characterized by their swirling nature,are mostly dominated by the evolution of vorticity. Small scale vortex structuresare essential to the look of smoke, fire, and related effects, whether producedby vortex interactions, jumps in density across interfaces (baroclinity), orviscous boundary layers. Classic Eulerian fluid solvers do not cost-effectively capturethese small-scale features in large domains. Lagrangian vortex methods showgreat promise from turbulence modelling, but face significant challenges in handlingboundary conditions, making them less attractive for computer graphics applications.This thesis proposes several novel solutions for the efficient simulationof small scale vortex features both in Eulerian and Lagrangian frameworks, extendingrobust Eulerian simulations with new algorithms inspired by Lagrangianvortex methods.
View record
Numerical solution of the time-harmonic Maxwell equations and incompressible magnetohydrodynamics problems (2011)
The goal of this thesis is to develop efficient numerical solvers for the time-harmonicMaxwell equations and for incompressible magnetohydrodynamics problems.The thesis consists of three components. In the first part, we present a fullyscalable parallel iterative solver for the time-harmonic Maxwell equations in mixedform with small wave numbers. We use the lowest order Nedelec elements of thefirst kind for the approximation of the vector field and standard nodal elements forthe Lagrange multiplier associated with the divergence constraint. The correspondinglinear system has a saddle point form, with inner iterations solved by preconditioned conjugate gradients. We demonstrate the performance ofour parallel solver on problems with constant and variable coefficients with up toapproximately 40 million degrees of freedom. Our numerical results indicate verygood scalability with the mesh size, on uniform, unstructured and locally refinedmeshes.In the second part, we introduce and analyze a mixed finite element method forthe numerical discretization of a stationary incompressible magnetohydrodynamicsproblem, in two and three dimensions. The velocity field is discretized usingdivergence-conforming Brezzi-Douglas-Marini (BDM) elements and the magneticfield is approximated by curl-conforming Nedelec elements. Key features of themethod are that it produces exactly divergence-free velocity approximations, andthat it correctly captures the strongest magnetic singularities in non-convex polyhedraldomains. We prove that the energy norm of the error is convergent in themesh size in general Lipschitz polyhedra under minimal regularity assumptions,and derive nearly optimal a-priori error estimates for the two-dimensional case. Wepresent a comprehensive set of numerical experiments, which indicate optimal convergence of the proposed method for two-dimensional as well as three-dimensionalproblems.Finally, in the third part we investigate preconditioned Krylov iterations forthe discretized stationary incompressible magnetohydrodynamics problems. Wepropose a preconditioner based on efficient preconditioners for the Maxwell andNavier-Stokes sub-systems. We show that many of the eigenvalues of the preconditionedsystem are tightly clustered, and hence, rapid convergence is accomplished.Our numerical results show that this approach performs quite well.
View record
Master's Student Supervision
Theses completed in 2010 or later are listed below. Please note that there is a 6-12 month delay to add the latest theses.
An experimental study of the incomplete Cholesky factorization for the anisotropic diffusion equation (2022)
We consider the numerical solution of large and sparse linear systems arising from a finite difference discretization of the anisotropic diffusion equation. We study preconditioned conjugate gradient as an iterative solver, with the incomplete Cholesky factorization as a preconditioner. The incomplete Cholesky factor is created using two different dynamic dropping parameters: a drop tolerance to control the magnitude of included values, and a fill factor to control the memory allocation/usage of the factor. For large matrices, we show that certain drop tolerances will achieve optimal results with a controllable amount of memory usage. Our numerical experiments are based on a fast C code that has been written as part of this research project, and which uses an asymptotically optimal version of the Cholesky factorization as its base. We illustrate our findings on linear systems of dimensions up to ten million degrees of freedom.
View record
Finite Difference Schemes for Elliptic Partial Differential Equations Requiring a Non-Uniform Mesh (2016)
A variety of finite difference schemes are explored for the numerical solution of elliptic partial differential equations, specifically the Poisson and convection-diffusion equations. Problems are investigated that require the use of a non-uniform or non-square mesh. This may be due to a non-square domain or a problem with a singularity. We explore the properties of the linear operators in the resulting systems of linear equations. In particular, we investigate the conditioning and eigenvalues of these operators, both numerically and in search of an approximation of these eigenvalues. We also investigate the choice of finite difference scheme with respect to accuracy and cost.
View record
Iterative Solution of a Mixed Finite Element Discretisation of an Incompressible Magnetohydrodynamics Problem (2014)
The aim of this thesis is to develop and numerically test a large scale preconditioned finite element implementation of an incompressible magnetohydrodynamics (MHD) model. To accomplish this, a broad-scope code has been generated using the finite element software package FEniCS and the linear algebra software PETSc. The code is modular, extremely flexible, and allows for implementing and testing different discretisations and linear algebra solvers with relatively modest effort. It can handle two- and three-dimensional problems in excess of 20 million degrees of freedom. Incompressible MHD describes the interaction between an incompressible electrically charged fluid governed by the incompressible Navier-Stokes equations coupled with electromagnetic effects from Maxwell’s equations in mixed form. We introduce a model problem and a mixed finite element discretisation based on using Taylor-Hood elements for the fluid variables and on a mixed N ́ed ́elec pair for the magnetic unknowns. We introduce three iteration strategies to handle the non-linearities present in the model, ranging from Picard iterations to completely decoupled schemes. Adapting and extending ideas introduced in [Dan Li, Numerical Solution of the Time-Harmonic Maxwell Equations and Incompressible Magnetohydrodynamics Problems, Ph.D. Dissertation, The University of British ii Columbia, 2010], we implement a preconditioning approach motivated by the block structure of the underlying linear systems in conjunction with state of the art preconditioners for the mixed Maxwell and Navier-Stokes subproblems. For the Picard iteration scheme we implement an inner-outer preconditioner. The numerical results presented in this thesis demonstrate the efficient performance of our preconditioned solution techniques and show good scalability with respect to the discretisation parameters.
View record
A parallel active-set method for solving frictional contact problems (2013)
Simulating frictional contact is a challenging computational task and there exist a variety of techniques to do so. One such technique, the staggered projections algorithm, requires the solution of two convex quadratic program (QP) subproblems at each iteration. We introduce a method, SCHURPA, which employs a primal-dual active-set strategy to efficiently solve these QPs based on a Schur-complement method. A single factorization of the initial saddle point system and a smaller dense Schur-complement is maintained to solve subsequent saddle point systems. Exploiting the parallelizability and warm-starting capabilities of the active-set method as well as the problem structure of the QPs yields a novel approach to the problem of frictional contact. Numerical results of a parallel GPU implementation using NVIDIA’s CUDA applied to a physical simulator of highly deformable bodies are presented.
View record
ILDL factorizations for sparse skew-symmetric matrices (2012)
Our goal is to solve a sparse skew-symmetric linear system efficiently. We propose a slight modification to the Bunch LDL\T factorization with partial pivoting strategy for skew-symmetric matrices, which saves approximately one third of the overall number of swaps of rows and columns on average. We also develop a rook pivoting strategy for this LDL\T factorization in Crout order. We derive an incomplete LDL\T factorization based on the full LDL\T factorization, with both partial pivoting strategy and rook pivoting strategy. The incomplete factorization can be used as a preconditioner for Krylov subspace solvers. Various column versions of ILUTP factorizations are investigated. We show that the approach taken saves approximately half of the work and storage compared to standard ILU factorizations that do not exploit the structure of the matrix, and the rook pivoting strategy is often superior to the partial pivoting strategy, in particular when the matrix is ill-conditioned and the overall number of elements in the rows of the factors is restricted.
View record
Spectral Properties of matrices Arising from Primal-Dual Interior-Point Methods for Convex Quadratic Programs (2012)
The solution of convex quadratic programs using primal-dual interior-point methods has at its core the solution of a series of linear systems, which are in practice commonly reduced by block Gaussian elimination from the original unsymmetric block 3-by-3 formulation to either a block 2-by-2 saddle-point matrix or a block 1-by-1 normal equations form. The 3-by-3 formulation can also be symmetrized with a similarity transformation if a symmetric solver is to be used. We examine whether this practice of reduction is beneficial for the solver. For each of these formulations we find analytical results for the following spectral properties: the inertia, condition number, conditions for nonsingularity, and bounds on the eigenvalues. While the reduced systems become increasingly ill-conditioned throughout the iterations except in special cases, the 3-by-3 formulations remain nonsingular and well-conditioned with only mild assumptions on the problem; with regularization the assumptions are further simplified. Numerical examples are used to support the analytical results. We conclude that the 3-by-3 formulations, unsymmetric or symmetric, unregularized or regularized, have superior spectral properties that support their use in practice.
View record
A complex analysis based derivation of multigrid smoothing factors of lexicographic Gauss-Seidel (2011)
This thesis aims to present a unified framework for deriving analytical formulas for multigrid smoothing factors in arbitrary dimensions, under certain simplifying assumptions. To derive these expressions we rely on complex analysis and geometric considerations, using the maximum modulus principle and Mobius transformations.We restrict our attention to pointwise and block lexicographic Gauss-Seidel smoothers on a d-dimensional uniform mesh, where the computational molecule of the associated discrete operator forms a 2d+1 point star. In the pointwise case the effect of a relaxation parameter, as well as different choices of mesh ratio, are analyzed. The results apply to any number of spatial dimensions, and are applicable to high-dimensional versions of a few common model problems with constant coefficients, including the Poisson and anisotropic diffusion equations as well as the convection-diffusion equation. We show that in most cases our formulas, exact under the simplifying assumptions of Local Fourier Analysis, form tight upper bounds for the asymptotic convergence of geometric multigrid in practice. We also show that there are asymmetric cases where lexicographic Gauss-Seidel smoothing outperforms red-black Gauss-Seidel smoothing; this occurs for certain model convection-diffusion equations with high mesh Reynolds numbers.
View record
Numerical Solution of Skew-Symmetric Linear Systems (2010)
We are concerned with iterative solvers for large and sparse skew-symmetric linear systems. First we discuss algorithms for computing incomplete factorizations as a source of preconditioners. This leads to a new Crout variant of Gaussian elimination for skew-symmetric matrices. Details on how to implement the algorithms efficiently are provided. A few numerical results are presented for these preconditioners. We also examine a specialized preconditioned minimum residual solver. An explicit derivation is given, detailing the effects of skew-symmetry on the algorithm.
View record
If this is your researcher profile you can log in to the Faculty & Staff portal to update your details and provide recruitment preferences.
Membership Status
Program Affiliations
Academic Unit(s)