Ivan G. Graham

Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens et al.

In this paper we prove, under mild conditions, that the positive definiteness of the circulant matrix appearing in the circulant embedding method is always guaranteed, provided the enclosing cube is sufficiently large. We examine in detail the case of the Matérn covariance, and prove (for fixed correlation length) that, as $h_0\rightarrow 0$, positive definiteness is guaranteed when the random field is sampled on a cube of size order $(1 + ν^{1/2} \log h_0^{-1})$ times larger than the size of the physical domain. (Here $h_0$ is the mesh spacing of the regular grid and $ν$ the Matérn smoothness parameter.) We show that the sampling cube can become smaller as the correlation length decreases when $h_0$ and $ν$ are fixed. Our results are confirmed by numerical experiments. We prove several results about the decay of the eigenvalues of the circulant matrix. These lead to the conjecture, verified by numerical experiment, that they decay with the same rate as the Karhunen--Loève eigenvalues of the covariance operator.

Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens et al.

In a previous paper (J. Comp. Phys. 230 (2011), 3668--3694), the authors proposed a new practical method for computing expected values of functionals of solutions for certain classes of elliptic partial differential equations with random coefficients. This method was based on combining quasi-Monte Carlo (QMC) methods for computing the expected values with circulant embedding methods for sampling the random field on a regular grid. It was found capable of handling fluid flow problems in random heterogeneous media with high stochastic dimension, but a convergence theory was missing. This paper provides a convergence analysis for the method in the case when the QMC method is a specially designed randomly shifted lattice rule. The convergence result depends on the eigenvalues of the underlying nested block circulant matrix and can be independent of the number of stochastic variables under certain assumptions. In fact the QMC analysis applies to general factorisations of the covariance matrix to sample the random field. The error analysis for the underlying fully discrete finite element method allows for locally refined meshes (via interpolation from a regular sampling grid of the random field). Numerical results on a non-regular domain with corner singularities in two spatial dimensions and on a regular domain in three spatial dimensions are included.

Timo Betcke, Simon N. Chandler-Wilde, Ivan G. Graham et al.

We consider the classical coupled, combined-field integral equation formulations for time-harmonic acoustic scattering by a sound soft bounded obstacle. In recent work, we have proved lower and upper bounds on the $L^2$ condition numbers for these formulations, and also on the norms of the classical acoustic single- and double-layer potential operators. These bounds to some extent make explicit the dependence of condition numbers on the wave number $k$, the geometry of the scatterer, and the coupling parameter. For example, with the usual choice of coupling parameter they show that, while the condition number grows like $k^{1/3}$ as $k\to\infty$, when the scatterer is a circle or sphere, it can grow as fast as $k^{7/5}$ for a class of `trapping' obstacles. In this paper we prove further bounds, sharpening and extending our previous results. In particular we show that there exist trapping obstacles for which the condition numbers grow as fast as $\exp(γk)$, for some $γ>0$, as $k\to\infty$ through some sequence. This result depends on exponential localisation bounds on Laplace eigenfunctions in an ellipse that we prove in the appendix. We also clarify the correct choice of coupling parameter in 2D for low $k$. In the second part of the paper we focus on the boundary element discretisation of these operators. We discuss the extent to which the bounds on the continuous operators are also satisfied by their discrete counterparts and, via numerical experiments, we provide supporting evidence for some of the theoretical results, both quantitative and asymptotic, indicating further which of the upper and lower bounds may be sharper.

Ivan G. Graham, Euan A. Spence, Eero Vainikko

In this paper we give new results on domain decomposition preconditioners for GMRES when computing piecewise-linear finite-element approximations of the Helmholtz equation $-Δu - (k^2+ {\rm i} \varepsilon)u = f$, with absorption parameter $\varepsilon \in \mathbb{R}$. Multigrid approximations of this equation with $\varepsilon \not= 0$ are commonly used as preconditioners for the pure Helmholtz case ($\varepsilon = 0$). However a rigorous theory for such (so-called "shifted Laplace") preconditioners, either for the pure Helmholtz equation, or even the absorptive equation ($\varepsilon \not=0$), is still missing. We present a new theory for the absorptive equation that provides rates of convergence for (left- or right-) preconditioned GMRES, via estimates of the norm and field of values of the preconditioned matrix. This theory uses a $k$- and $\varepsilon$-explicit coercivity result for the underlying sesquilinear form and shows, for example, that if $|\varepsilon|\sim k^2$, then classical overlapping additive Schwarz will perform optimally for the absorptive problem, provided the subdomain and coarse mesh diameters are carefully chosen. Extensive numerical experiments are given that support the theoretical results. The theory for the absorptive case gives insight into how its domain decomposition approximations perform as preconditioners for the pure Helmholtz case $\varepsilon = 0$. At the end of the paper we propose a (scalable) multilevel preconditioner for the pure Helmholtz problem that has an empirical computation time complexity of about $\mathcal{O}(n^{4/3})$ for solving finite element systems of size $n=\mathcal{O}(k^3)$, where we have chosen the mesh diameter $h \sim k^{-3/2}$ to avoid the pollution effect. Experiments on problems with $h\sim k^{-1}$, i.e. a fixed number of grid points per wavelength, are also given.

Marcella Bonazzoli, Victorita Dolean, Ivan G. Graham et al.

This paper rigorously analyses preconditioners for the time-harmonic Maxwell equations with absorption, where the PDE is discretised using curl-conforming finite-element methods of fixed, arbitrary order and the preconditioner is constructed using Additive Schwarz domain decomposition methods. The theory developed here shows that if the absorption is large enough, and if the subdomain and coarse mesh diameters and overlap are chosen appropriately, then the classical two-level overlapping Additive Schwarz preconditioner (with PEC boundary conditions on the subdomains) performs optimally -- in the sense that GMRES converges in a wavenumber-independent number of iterations -- for the problem with absorption. An important feature of the theory is that it allows the coarse space to be built from low-order elements even if the PDE is discretised using high-order elements. It also shows that additive methods with minimal overlap can be robust. Numerical experiments are given that illustrate the theory and its dependence on various parameters. These experiments motivate some extensions of the preconditioners which have better robustness for problems with less absorption, including the propagative case. At the end of the paper we illustrate the performance of these on two substantial applications; the first (a problem with absorption arising from medical imaging) shows the empirical robustness of the preconditioner against heterogeneity, and the second (scattering by a COBRA cavity) shows good scalability of the preconditioner with up to 3,000 processors.

Théophile Chaumont-Frelet, Victorita Dolean, Mark Fry et al.

GenEO (`Generalised Eigenvalue problems on the Overlap') is a method for constructing coarse spaces used in the preconditioning of iterative solvers for discrete PDEs. This method combines a (small) number of modes of local PDE eigenproblems to obtain a global coarse space. A coarse solve is then combined with local solves of the global PDE to obtain the preconditioner. A substantial theory for GenEO has been developed for the case when the local elgenproblems are positive semi-definite. This has been applied mostly to positive definite global PDEs, but also recently extended to the case of convection--diffusion--reaction problems, which may be neither self-adjoint, nor positive definite. However, when the global problem is highly indefinite, coarse spaces built from positive semi-definite local eigenproblems fail to be robust in practice. In this paper we consider highly indefinite global PDE problems, characterised by a large parameter $k$ (allowing also highly variable coefficients), and we develop a new spectral coarse space built from solving eigenvalue problems based on \textit{local copies of the global problem}. We put no constraint on the diameters of the local domains, thus allowing the local eigenvalue problems to be indefinite. The new method (which we call $H_k$-GenEO) is seen to be much more robust as $k$ increases than methods based on positive semi-definite eigenproblems. We provide sufficient conditions for robustness of the preconditioned GMRES iterative method, in terms of the tolerance of the local eigenproblems and the size of the subdomains for the local PDE solves. In practice the method is observed to be robust with respect to $k$ under even weaker conditions on the local eigenproblem tolerance. The experiments also suggest the method can be resilient to high variation in PDE coefficients.

Ivan G. Graham, Matthew J. Parkinson, Robert Scheichl

We describe modern variants of Monte Carlo methods for Uncertainty Quantification (UQ) of the Neutron Transport Equation, when it is approximated by the discrete ordinates method with diamond differencing. We focus on the mono-energetic 1D slab geometry problem, with isotropic scattering, where the cross-sections are log-normal correlated random fields of possibly low regularity. The paper includes an outline of novel theoretical results on the convergence of the discrete scheme, in the cases of both spatially variable and random cross-sections. We also describe the theory and practice of algorithms for quantifying the uncertainty of a linear functional of the scalar flux, using Monte Carlo and quasi-Monte Carlo methods, and their multilevel variants. A hybrid iterative/direct solver for computing each realisation of the functional is also presented. Numerical experiments show the effectiveness of the hybrid solver and the gains that are possible through quasi-Monte Carlo sampling and multilevel variance reduction. For the multilevel quasi-Monte Carlo method, we observe gains in the computational $\varepsilon$-cost of up to 2 orders of magnitude over the standard Monte Carlo method, and we explain this theoretically. Experiments on problems with up to several thousand stochastic dimensions are included.

Alexander D. Gilbert, Ivan G. Graham, Frances Y. Kuo et al.

We consider the forward problem of uncertainty quantification for the generalised Dirichlet eigenvalue problem for a coercive second order partial differential operator with random coefficients, motivated by problems in structural mechanics, photonic crystals and neutron diffusion. The PDE coefficients are assumed to be uniformly bounded random fields, represented as infinite series parametrised by uniformly distributed i.i.d. random variables. The expectation of the fundamental eigenvalue of this problem is computed by (a) truncating the infinite series which define the coefficients; (b) approximating the resulting truncated problem using lowest order conforming finite elements and a sparse matrix eigenvalue solver; and (c) approximating the resulting finite (but high dimensional) integral by a randomly shifted quasi-Monte Carlo lattice rule, with specially chosen generating vector. We prove error estimates for the combined error, which depend on the truncation dimension $s$, the finite element mesh diameter $h$, and the number of quasi-Monte Carlo samples $N$. Under suitable regularity assumptions, our bounds are of the particular form $\mathcal{O}(h^2+N^{-1+δ})$, where $δ>0$ is arbitrary and the hidden constant is independent of the truncation dimension, which needs to grow as $h\to 0$ and $N\to\infty$. Although the eigenvalue problem is nonlinear, which means it is generally considered harder than the analogous source problem, in almost all cases we obtain error bounds that converge at the same rate as the corresponding rate for the source problem. The proof involves a detailed study of the regularity of the fundamental eigenvalue as a function of the random parameters. As a key intermediate result in the analysis, we prove that the spectral gap (between the fundamental and the second eigenvalues) is uniformly positive over all realisations of the random problem.