Past Computational Mathematics and Applications Seminar

Following the framework of the heterogeneous multiscale method, we present a numerical method for nonlinear elliptic homogenization problems. We briefly review the numerical, relying on an efficient coupling of macro and micro solvers, for linear problems. A fully discrete analysis is then given for nonlinear (nonmonotone) problems, optimal convergence rates in the H¹ and L² norms are derived and the uniqueness of the method is shown on sufficiently fine macro and micro meshes.

Numerical examples confirm the theoretical convergence rates and illustrate the performance and versatility of our approach.

This talk is about the Induced Dimension Reduction (IDR) methods developed by Peter Sonneveld and, more recently, Martin van Gijzen. We sketch the history, outline the underlying principle, and give a few details about different points of view on this class of Krylov subspace methods. If time permits, we briefly outline some recent developments in this field and the benefits and drawbacks of these and IDR methods in general.

Many continuum models have been derived in recent years which describe the self-assembly of industrially utilisable crystalline films to a level of detail that allows qualitative comparisons with experiments. For thin-film problems, where the characteristic length scales in vertical and horizontal directions differ significantly, the governing surface diffusion equations can be reduced to simpler PDEs by making use of asymptotic expansions. Many mathematical problems and solutions emerge from such new evolution equations and many of them remind of Cahn-Hilliard type equations. The surface diffusion models are of high, of fourth or even sixth, order.

We present the modeling, model reduction and simulation results for heteroepitaxial growth as for Ge/Si quantum dot self-assembly. The numerical methods we are using are based on trigonometric interpolation. These kind of pseudospectral methods seem very well suited for simulating the coarsening of large quantum dot arrays. When the anisotropy of the growing crystalline film is strong, it might become necessary to add a corner regularisation to the model. Then the transition region between neighboring facets is still smooth, but its scale is rather small. In this case it might be useful to think about an adaptive extension of the existing method.

Figure 1: Ostwald ripening process of quantum dots depicted at consecutive time points. One fourth of the whole, periodic, simulated domain is shown.

Joint work with Peter Evans and Barbara Wagner

This seminar will be held at the Rutherford Appleton Laboratory near Didcot.

Abstract:

Numerical calculations of laminar flow in a two-dimensional channel with a sudden expansion exhibit a symmetry-breaking bifurcation at Reynolds number 40.45 when the expansion ratio is 3:1. In the experiments reported by Fearn, Mullin and Cliffe [1] there is a large perturbation to this bifurcation and the agreement with the numerical calculations is surprisingly poor. Possible reasons for this discrepancy are explored using modern techniques for uncertainty quantification.

When experimental equipment is constructed there are, inevitably, small manufacturing imperfections that can break the symmetry in the apparatus. In this work we considered a simple model for these imperfections. It was assumed that the inlet section of the channel was displaced by a small amount and that the centre line of the inlet section was not parallel to the centre line of the outlet section. Both imperfections were modelled as normal random variables with variance equal to the manufacturing tolerance. Thus the problem to be solved is the Navier-Stokes equations in a geometry with small random perturbations. A co-ordinate transformation technique was used to transform the problem to a fixed deterministic domain but with random coefficient appearing in the transformed Navier-Stokes equations. The resulting equations were solved using a stochastic collocation technique that took into account the fact that the problem has a discontinuity in parameter space arising from the bifurcation structure in the problem.

The numerical results are in the form of an approximation to a probability measure on the set of bifurcation diagrams. The experimental data of Fearn, Mullin and Cliffe are consistent with the computed solutions, so it appears that a satisfactory explanation for the large perturbation can be provided by manufacturing imperfections in the experimental apparatus.

The work demonstrates that modern methods for uncertainty quantification can be applied successfully to a bifurcation problem arising in fluid mechanics. It should be possible to apply similar techniques to a wide range of bifurcation problems in fluid mechanics in the future.

References:

[1] R M Fearn, T Mullin and K A Cliffe Nonlinear flow phenomena in a symmetric sudden expansion, J. Fluid Mech. 211, 595-608, 1990.

Please note that this is a short notice change from the originally advertised talk by Dr Shahrokh Shahpar (Rolls-Royce plc.)

The new talk "Multilevel Monte Carlo method" is given by Mike Giles, Oxford-Man Institute of Quantitative Finance, Mathematical Institute, University of Oxford.

Joint work with Rob Scheichl, Aretha Teckentrup (Bath) and Andrew Cliffe (Nottingham)

Over the last several decades, many mesh generation methods and a plethora of adaptive methods for solving differential equations have been developed.  In this talk, we take a general approach for describing the mesh generation problem, which can be considered as being in some sense equivalent to determining a coordinate transformation between physical space and a computational space.  Our description provides some new theoretical insights into precisely what is accomplished from mesh equidistribution (which is a standard adaptivity tool used in practice) and mesh alignment.  We show how variational mesh generation algorithms, which have historically been the most common and important ones, can generally be compared using these mesh generation principles.  Lastly, we relate these to a variety of moving mesh methods for solving time-dependent PDEs.

This is joint work with Weizhang Huang, Kansas University

In this talk we discuss the design of efficient numerical methods for solving symmetric indefinite linear systems arising from mixed approximation of elliptic PDE problems with associated constraints. Examples include linear elasticity (Navier-Lame equations), steady fluid flow (Stokes' equations) and electromagnetism (Maxwell's equations).

The novel feature of our iterative solution approach is the incorporation of error control in the natural "energy" norm in combination with an a posteriori estimator for the PDE approximation error. This leads to a robust and optimally efficient stopping criterion: the iteration is terminated as soon as the algebraic error is insignificant compared to the approximation error. We describe a "proof of concept" MATLAB implementation of this algorithm, which we call EST_MINRES, and we illustrate its effectiveness when integrated into our Incompressible Flow Iterative Solution Software (IFISS) package (http://www.manchester.ac.uk/ifiss/).

In this talk, we focus on the analytical properties of a decentralized auction, namely the PAUSE Auction Procedure. We prove that the revenue of the auctioneer from PAUSE is greater than or equal to the profit from the well-known VCG auction when there are only two bidders and provide lower bounds on the profit for arbitrary number of bidders. Based on these bounds and observations from auctions with few items, we propose a modification of the procedure that increases the profit. We believe that this study, which is still in progress, will be a milestone in designing better decentralized auctions since it is the first analytical study on such auctions with promising results.

Polynomial Programs are ussually solved by using hierarchies of convex relaxations. This scheme rapidly becomes computationally expensive and is often tractable only for problems of small sizes. We propose an iterative scheme that improves an initial relaxation without incurring exponential growth in size. The key ingredient is a dynamic scheme for generating valid polynomial inequalities for general polynomial programs. These valid inequalities are then used to construct better approximations of the original problem. As a result, the proposed scheme is in principle scalable to large general combinatorial optimization problems.

Joint work with Bissan Ghaddar and Miguel Anjos

Based on an MPI library written over 10 years ago, OP2 is a new open-source library which is aimed at application developers using unstructured grids. Using a single API, it targets a variety of HPC architectures, including both manycore GPUs and multicore CPUs with vector units. The talk will cover the API design, key aspects of the parallel implementation on the different platforms, and preliminary performance results on a small but representative CFD test code.

Project homepage: http://people.maths.ox.ac.uk/gilesm/op2/

We consider the iterative solution of large sparse linear least squares (LS) problems. Specifically, we focus on the design and implementation of reliable stopping criteria for the widely-used algorithm LSQR of Paige and Saunders. First we perform a backward perturbation analysis of the LS problem. We show why certain projections of the residual vector are good measures of convergence, and we propose stopping criteria that use these quantities. These projections are too expensive to compute to be used directly in practice. We show how to estimate them efficiently at every iteration of the algorithm LSQR. Our proposed stopping criteria can therefore be used in practice.

This talk is based on joint work with Xiao-Wen Chang, Chris Paige, Pavel Jiranek, and Serge Gratton.

In various fields of application, one must solve very large scale boundary value problems using parallel solvers and supercomputers. The domain decomposition approach partitions the large computational domain into smaller computational subdomains. In order to speed up the convergence, we have several ``optimized'' algorithm that use Robin transmission conditions across the artificial interfaces (FETI-2LM). It is known that this approach alone is not sufficient: as the number of subdomains increases, the number of iterations required for convergence also increases and hence the parallel speedup is lost. A known solution for classical Schwarz methods as well as FETI algorithms is to incorporate a ``coarse grid correction'', which is able to transmit low-frequency information more quickly across the whole domain. Such algorithms are known to ``scale weakly'' to large supercomputers. A coarse grid correction is also necessary for FETI-2LM methods. In this talk, we will introduce and analyze coarse grid correction algorithms for FETI-2LM domain decomposition methods.

Can one of the most important Italian Renaissance frescoes reduced in hundreds of thousand fragments by a bombing during the Second World War be re-composed after more than 60 years from its damage? Can we reconstruct the missing parts and can we say something about their original color?

In this talk we would like to exemplify, hopefully effectively by taking advantage of the seduction of art, how mathematics today can be applied in real-life problems which were considered unsolvable only few years ago.

Implementations of the revised simplex method for solving large scale sparse linear programming (LP) problems are highly efficient for single-core architectures. This talk will discuss the limitations of the underlying techniques in the context of modern multi-core architectures, in particular with respect to memory access. Novel techniques for implementing the dual revised simplex method will be introduced, and their use in developing a dual revised simplex solver for multi-core architectures will be described.

We propose and analyze a primal-dual active set method for local and non-local vector-valued Allen-Cahn variational inequalities.

We show existence and uniqueness of a solution for the non-local vector-valued Allen-Cahn variational inequality in a formulation involving Lagrange multipliers for local and non-local constraints. Furthermore, convergence of the algorithm is shown by interpreting the approach as a semi-smooth Newton method and numerical simulations are presented.

Optimization problems with time-periodic parabolic PDE constraints can arise in important chemical engineering applications, e.g., in periodic adsorption processes. I will present a novel direct numerical method for this problem class. The main numerical challenges are the high nonlinearity and high dimensionality of the discretized problem. The method is based on Direct Multiple Shooting and inexact Sequential Quadratic Programming with globalization of convergence based on natural level functions. I will highlight the use of a generalized Richardson iteration with a novel two-grid Newton-Picard preconditioner for the solution of the quadratic subproblems. At the end of the talk I will explain the principle of Simulated Moving Bed processes and conclude with numerical results for optimization of such a process.

Efficient linear and infinitely smooth approximation of functions from equidistant samples is a fascinating problem, at least since Runge showed in 1901 that it is not delivered by the interpolating polynomial.

In 1988, I suggested to substitute linear rational for polynomial interpolation by replacing the denominator 1 with a polynomial depending on the nodes, though not on the interpolated function. Unfortunately the so-obtained interpolant converges merely as the square of the mesh size. In 2007, Floater and Hormann have given for every integer a denominator that yields convergence of that prescribed order.

In the present talk I shall present the corresponding interpolant as well as some of its applications to differentiation, integration and the solution of boundary value problems. This is joint work with Georges Klein and Michael Floater.

The Bi-Conjugate Gradient method (BiCG) is a well-known iterative solver (Krylov method) for linear systems of equations, proposed about 35 years ago, and the basis for some of the most successful iterative methods today, like BiCGSTAB. Nevertheless, the convergence behavior is poorly understood. The method satisfies a Petrov-Galerkin property, and hence its residual is constrained to a space of decreasing dimension (decreasing one per iteration). However, that does not explain why, for many problems, the method converges in, say, a hundred or a few hundred iterations for problems involving a hundred thousand or a million unknowns. For many problems, BiCG converges not much slower than an optimal method, like GMRES, even though the method does not satisfy any optimality properties. In fact, Anne Greenbaum showed that every three-term recurrence, for the first (n/2)+1 iterations (for a system of dimension n), is BiCG for some initial 'left' starting vector. So, why does the method work so well in most cases? We will introduce Krylov methods, discuss the convergence of optimal methods, describe the BiCG method, and provide an analysis of its convergence behavior.

Algebraic multigrid methods are nowadays popular to solve linear systems arising from the discretization of elliptic PDEs. They try to combine the efficiency of well tuned specific schemes like classical (geometric-based) multigrid methods, with the ease of use of general purpose preconditioning techniques. This requires to define automatic coarsening procedures, which set up an hierarchy of coarse representations of the problem at hand using only information from the system matrix.

In this talk, we focus on aggregation-based algebraic multigrid methods. With these, the coarse unknowns are simply formed by grouping variables in disjoint subset called aggregates.

In the first part of the talk, we consider symmetric M-matrices with nonnegative row-sum. We show how aggregates can then be formed in such a way that the resulting method satisfies a prescribed bound on its convergence rate. That is, instead of the classical paradigm that applies an algorithm and then performs its analysis, the analysis is integrated in the set up phase so as to enforce minimal quality requirements. As a result, we obtain one of the first algebraic multigrid method with full convergence proof. The efficiency of the method is further illustrated by numerical results performed on finite difference or linear finite element discretizations of second order elliptic PDEs; the set of problems includes problems with jumps, anisotropy, reentering corners and/or unstructured meshes, sometimes with local refinement.

In the second part of the talk, we discuss nonsymmetric problems. We show how the previous approach can be extended to M-matrices with row- and column-sum both nonnegative, which includes some stable discretizations of convection-diffusion equations with divergence free convective flow. Some (preliminary) numerical results are also presented.

This is joint work with Artem Napov.

Starting from a Navier-Stokes-Cahn-Hilliard equation for a two-phase flow problem we discuss efficient numerical approaches based on adaptive finite element methods. Various extensions of the model are discussed: a) we consider the model on implicitly described geometries, which is used to simulate the sliding of droplets over nano-patterned surfaces, b) we consider the effect of soluble surfactants and show its influence on tip splitting of droplets under shear flow, and c) we consider bijels as a new class of soft matter materials, in which colloidal particles are jammed on the fluid-fluid interface and effect the motion of the interface due to an elastic force.

The work is based on joint work with Sebastian Aland (TU Dresden), John Lowengrub (UC Irvine) and Knut Erik Teigen (U Trondheim).

We extend for the first time the linear discretization theory of Schaback, developed for meshfree methods, to nonlinear operator equations, relying heavily on methods of Böhmer, Vol I. There is no restriction to elliptic problems or to symmetric numerical methods like Galerkin techniques.

Trial spaces can be arbitrary, but have to approximate the solution well, and testing can be weak or strong. We present Galerkin techniques as an example. On the downside, stability is not easy to prove for special applications, and numerical methods have to be formulated as optimization problems. Results of this discretization theory cover error bounds and convergence rates. These results remain valid for the general case of fully nonlinear elliptic differential equations of second order. Some numerical examples are added for illustration.

We present an O(N logN) algorithm for the calculation of the first N coefficients in an expansion of an analytic function in Legendre polynomials. In essence, the algorithm consists of an integration of a suitably weighted function along an ellipse, a task which can be accomplished with Fast Fourier Transform, followed by some post-processing.