Showing 1–50 of 62 results for author: Scheichl, R

A Mixed Multiscale Spectral Generalized Finite Element Method

Authors: Christian Alber, Chupeng Ma, Robert Scheichl

Abstract: We present a multiscale mixed finite element method for solving second order elliptic equations with general $L^{\infty}$-coefficients arising from flow in highly heterogeneous porous media. Our approach is based on a multiscale spectral generalized finite element method (MS-GFEM) and exploits the superior local mass conservation properties of mixed finite elements. Following the MS-GFEM framework… ▽ More We present a multiscale mixed finite element method for solving second order elliptic equations with general $L^{\infty}$-coefficients arising from flow in highly heterogeneous porous media. Our approach is based on a multiscale spectral generalized finite element method (MS-GFEM) and exploits the superior local mass conservation properties of mixed finite elements. Following the MS-GFEM framework, optimal local approximation spaces are built for the velocity field by solving local eigenvalue problems over generalized harmonic spaces. The resulting global velocity space is then enriched suitably to ensure inf-sup stability. We develop the mixed MS-GFEM for both continuous and discrete formulations, with Raviart-Thomas based mixed finite elements underlying the discrete method. Exponential convergence with respect to local degrees of freedom is proven at both the continuous and discrete levels. Numerical results are presented to support the theory and to validate the proposed method. △ Less

Submitted 4 April, 2024; v1 submitted 25 March, 2024; originally announced March 2024.

Democratizing Uncertainty Quantification

Authors: Linus Seelinger, Anne Reinarz, Mikkel B. Lykkegaard, Robert Akers, Amal M. A. Alghamdi, David Aristoff, Wolfgang Bangerth, Jean Bénézech, Matteo Diez, Kurt Frey, John D. Jakeman, Jakob S. Jørgensen, Ki-Tae Kim, Massimiliano Martinelli, Matthew Parno, Riccardo Pellegrini, Noemi Petra, Nicolai A. B. Riis, Katherine Rosenfeld, Andrea Serani, Lorenzo Tamellini, Umberto Villa, Tim J. Dodwell, Robert Scheichl

Abstract: Uncertainty Quantification (UQ) is vital to safety-critical model-based analyses, but the widespread adoption of sophisticated UQ methods is limited by technical complexity. In this paper, we introduce UM-Bridge (the UQ and Modeling Bridge), a high-level abstraction and software protocol that facilitates universal interoperability of UQ software with simulation codes. It breaks down the technical… ▽ More Uncertainty Quantification (UQ) is vital to safety-critical model-based analyses, but the widespread adoption of sophisticated UQ methods is limited by technical complexity. In this paper, we introduce UM-Bridge (the UQ and Modeling Bridge), a high-level abstraction and software protocol that facilitates universal interoperability of UQ software with simulation codes. It breaks down the technical complexity of advanced UQ applications and enables separation of concerns between experts. UM-Bridge democratizes UQ by allowing effective interdisciplinary collaboration, accelerating the development of advanced UQ methods, and making it easy to perform UQ analyses from prototype to High Performance Computing (HPC) scale. In addition, we present a library of ready-to-run UQ benchmark problems, all easily accessible through UM-Bridge. These benchmarks support UQ methodology research, enabling reproducible performance comparisons. We demonstrate UM-Bridge with several scientific applications, harnessing HPC resources even using UQ codes not designed with HPC support. △ Less

Submitted 9 May, 2024; v1 submitted 21 February, 2024; originally announced February 2024.

Comments: Fix broken reference, add R. Akers as author for contribution to tritium benchmark

Tractable Optimal Experimental Design using Transport Maps

Authors: Karina Koval, Roland Herzog, Robert Scheichl

Abstract: We present a flexible method for computing Bayesian optimal experimental designs (BOEDs) for inverse problems with intractable posteriors. The approach is applicable to a wide range of BOED problems and can accommodate various optimality criteria, prior distributions and noise models. The key to our approach is the construction of a transport-map-based surrogate to the joint probability law of the… ▽ More We present a flexible method for computing Bayesian optimal experimental designs (BOEDs) for inverse problems with intractable posteriors. The approach is applicable to a wide range of BOED problems and can accommodate various optimality criteria, prior distributions and noise models. The key to our approach is the construction of a transport-map-based surrogate to the joint probability law of the design, observational and inference random variables. This order-preserving transport map is constructed using tensor trains and can be used to efficiently sample from (and evaluate approximate densities of) conditional distributions that are required in the evaluation of many commonly-used optimality criteria. The algorithm is also extended to sequential data acquisition problems, where experiments can be performed in sequence to update the state of knowledge about the unknown parameters. The sequential BOED problem is made computationally feasible by preconditioning the approximation of the joint density at the current stage using transport maps constructed at previous stages. The flexibility of our approach in finding optimal designs is illustrated with some numerical examples inspired by disease modeling and the reconstruction of subsurface structures in aquifers. △ Less

Submitted 22 April, 2024; v1 submitted 15 January, 2024; originally announced January 2024.

MSC Class: 62K05; 62F15; 65K10; 65L09; 65N21; 15A69

A complex-projected Rayleigh quotient iteration for targeting interior eigenvalues

Authors: Nils Friess, Alexander D. Gilbert, Robert Scheichl

Abstract: We introduce a new Projected Rayleigh Quotient Iteration aimed at improving the convergence behaviour of classic Rayleigh Quotient iteration (RQI) by incorporating approximate information about the target eigenvector at each step. While classic RQI exhibits local cubic convergence for Hermitian matrices, its global behaviour can be unpredictable, whereby it may converge to an eigenvalue far away f… ▽ More We introduce a new Projected Rayleigh Quotient Iteration aimed at improving the convergence behaviour of classic Rayleigh Quotient iteration (RQI) by incorporating approximate information about the target eigenvector at each step. While classic RQI exhibits local cubic convergence for Hermitian matrices, its global behaviour can be unpredictable, whereby it may converge to an eigenvalue far away from the target, even when started with accurate initial conditions. This problem is exacerbated when the eigenvalues are closely spaced. The key idea of the new algorithm is at each step to add a complex-valued projection to the original matrix (that depends on the current eigenvector approximation), such that the unwanted eigenvalues are lifted into the complex plane while the target stays close to the real line, thereby increasing the spacing between the target eigenvalue and the rest of the spectrum. Making better use of the eigenvector approximation leads to more robust convergence behaviour and the new method converges reliably to the correct target eigenpair for a significantly wider range of initial vectors than does classic RQI. We prove that the method converges locally cubically and we present several numerical examples demonstrating the improved global convergence behaviour. In particular, we apply it to compute eigenvalues in a band-gap spectrum of a Sturm-Liouville operator used to model photonic crystal fibres, where the target and unwanted eigenvalues are closely spaced. The examples show that the new method converges to the desired eigenpair even when the eigenvalue spacing is very small, often succeeding when classic RQI fails. △ Less

Submitted 6 December, 2023; v1 submitted 5 December, 2023; originally announced December 2023.

Comments: 23 pages, 6 figures, 1 table

Improving Met Office Weather and Climate Forecasts with Bespoke Multigrid Solvers

Authors: Andrew Malcolm, Eike H. Müller, Robert Scheichl

Abstract: At the heart of the Met Office climate and weather forecasting capabilities lies a sophisticated numerical model which solves the equations of large-scale atmospheric flow. Since this model uses semi-implicit time-stepping, it requires the repeated solution of a large sparse system of linear equations with hundreds of millions of unknowns. This is one of the computational bottlenecks of operationa… ▽ More At the heart of the Met Office climate and weather forecasting capabilities lies a sophisticated numerical model which solves the equations of large-scale atmospheric flow. Since this model uses semi-implicit time-stepping, it requires the repeated solution of a large sparse system of linear equations with hundreds of millions of unknowns. This is one of the computational bottlenecks of operational forecasts and efficient numerical algorithms are crucial to ensure optimal performance. We developed and implemented a bespoke multigrid solver to address this challenge. Our solver reduces the time for solving the linear system by a factor two, compared to the previously used BiCGStab method. This leads to significant improvements of overall model performance: global forecasts can be produced 10%-15% faster. Multigrid also avoids stagnating convergence of the iterative scheme in single precision. By allowing better utilisation of computational resources, our work has led to estimated annual cost savings of GBP 300k for the Met Office. △ Less

Submitted 23 November, 2023; v1 submitted 10 July, 2023; originally announced July 2023.

Comments: 6 pages, 2 figures, 1 table; This is a preprint of the following chapter: A. Malcolm, E. H. Mueller, R. Scheichl: "Improving Met Office Weather and Climate Forecasts with Bespoke Multigrid Solvers", to be published in "More UK Success Stories in Industrial Mathematics", edited by Philip J. Aston, 2024, Springer reproduced with permission of Springer Nature Switzerland AG

MSC Class: 65F08; 65N55; 65Y05; 76M20; 86A10 ACM Class: G.1.3; G.1.8; J.2

Lowering the Entry Bar to HPC-Scale Uncertainty Quantification

Authors: Linus Seelinger, Anne Reinarz, Jean Benezech, Mikkel Bue Lykkegaard, Lorenzo Tamellini, Robert Scheichl

Abstract: Treating uncertainties in models is essential in many fields of science and engineering. Uncertainty quantification (UQ) on complex and computationally costly numerical models necessitates a combination of efficient model solvers, advanced UQ methods and HPC-scale resources. The resulting technical complexities as well as lack of separation of concerns between UQ and model experts is holding back… ▽ More Treating uncertainties in models is essential in many fields of science and engineering. Uncertainty quantification (UQ) on complex and computationally costly numerical models necessitates a combination of efficient model solvers, advanced UQ methods and HPC-scale resources. The resulting technical complexities as well as lack of separation of concerns between UQ and model experts is holding back many interesting UQ applications. The aim of this paper is to close the gap between advanced UQ methods and advanced models by removing the hurdle of complex software stack integration, which in turn will offer a straightforward way to scale even prototype-grade UQ applications to high-performance resources. We achieve this goal by introducing a parallel software architecture based on UM-Bridge, a universal interface for linking UQ and models. We present three realistic applications from different areas of science and engineering, scaling from single machines to large clusters on the Google Cloud Platform. △ Less

Submitted 27 April, 2023; originally announced April 2023.

Multilevel Monte Carlo methods for stochastic convection-diffusion eigenvalue problems

Authors: Tiangang Cui, Hans De Sterck, Alexander D. Gilbert, Stanislav Polishchuk, Robert Scheichl

Abstract: We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection-diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh q… ▽ More We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection-diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov--Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases. △ Less

Submitted 12 February, 2024; v1 submitted 7 March, 2023; originally announced March 2023.

MSC Class: 65N25

Scalable multiscale-spectral GFEM with an application to composite aero-structures

Authors: Jean Bénézech, Linus Seelinger, Peter Bastian, Richard Butler, Timothy Dodwell, Chupeng Ma, Robert Scheichl

Abstract: In this paper, the first large-scale application of multiscale-spectral generalized finite element methods (MS-GFEM) to composite aero-structures is presented. The crucial novelty lies in the introduction of A-harmonicity in the local approximation spaces, which in contrast to [Babuska, Lipton, Multiscale Model. Simul. 9, 2011] is enforced more efficiently via a constraint in the local eigenproble… ▽ More In this paper, the first large-scale application of multiscale-spectral generalized finite element methods (MS-GFEM) to composite aero-structures is presented. The crucial novelty lies in the introduction of A-harmonicity in the local approximation spaces, which in contrast to [Babuska, Lipton, Multiscale Model. Simul. 9, 2011] is enforced more efficiently via a constraint in the local eigenproblems. This significant modification leads to excellent approximation properties, which turn out to be essential to capture accurately material strains and stresses with a low dimensional approximation space, hence maximising model order reduction. The implementation of the framework in the DUNE software package, as well as a detailed description of all components of the method are presented and exemplified on a composite laminated beam under compressive loading. The excellent parallel scalability of the method, as well as its superior performance compared to the related, previously introduced GenEO method are demonstrated on two realistic application cases, including a C-shaped wing spar with complex geometry. Further, by allowing low-cost approximate solves for closely related models or geometries this efficient, novel technology provides the basis for future applications in optimisation or uncertainty quantification on challenging problems in composite aero-structures. △ Less

Submitted 1 March, 2023; v1 submitted 24 November, 2022; originally announced November 2022.

Deep importance sampling using tensor trains with application to a priori and a posteriori rare event estimation

Authors: Tiangang Cui, Sergey Dolgov, Robert Scheichl

Abstract: We propose a deep importance sampling method that is suitable for estimating rare event probabilities in high-dimensional problems. We approximate the optimal importance distribution in a general importance sampling problem as the pushforward of a reference distribution under a composition of order-preserving transformations, in which each transformation is formed by a squared tensor-train decompo… ▽ More We propose a deep importance sampling method that is suitable for estimating rare event probabilities in high-dimensional problems. We approximate the optimal importance distribution in a general importance sampling problem as the pushforward of a reference distribution under a composition of order-preserving transformations, in which each transformation is formed by a squared tensor-train decomposition. The squared tensor-train decomposition provides a scalable ansatz for building order-preserving high-dimensional transformations via density approximations. The use of composition of maps moving along a sequence of bridging densities alleviates the difficulty of directly approximating concentrated density functions. To compute expectations over unnormalized probability distributions, we design a ratio estimator that estimates the normalizing constant using a separate importance distribution, again constructed via a composition of transformations in tensor-train format. This offers better theoretical variance reduction compared with self-normalized importance sampling, and thus opens the door to efficient computation of rare event probabilities in Bayesian inference problems. Numerical experiments on problems constrained by differential equations show little to no increase in the computational complexity with the event probability going to zero, and allow to compute hitherto unattainable estimates of rare event probabilities for complex, high-dimensional posterior densities. △ Less

Submitted 24 May, 2023; v1 submitted 5 September, 2022; originally announced September 2022.

Adaptive multilevel subset simulation with selective refinement

Authors: Daniel Elfverson, Robert Scheichl, Simon Weissmann, F. Alejandro DiazDelaO

Abstract: In this work we propose an adaptive multilevel version of subset simulation to estimate the probability of rare events for complex physical systems. Given a sequence of nested failure domains of increasing size, the rare event probability is expressed as a product of conditional probabilities. The proposed new estimator uses different model resolutions and varying numbers of samples across the hie… ▽ More In this work we propose an adaptive multilevel version of subset simulation to estimate the probability of rare events for complex physical systems. Given a sequence of nested failure domains of increasing size, the rare event probability is expressed as a product of conditional probabilities. The proposed new estimator uses different model resolutions and varying numbers of samples across the hierarchy of nested failure sets. In order to dramatically reduce the computational cost, we construct the intermediate failure sets such that only a small number of expensive high-resolution model evaluations are needed, whilst the majority of samples can be taken from inexpensive low-resolution simulations. A key idea in our new estimator is the use of a posteriori error estimators combined with a selective mesh refinement strategy to guarantee the critical subset property that may be violated when changing model resolution from one failure set to the next. The efficiency gains and the statistical properties of the estimator are investigated both theoretically via shaking transformations, as well as numerically. On a model problem from subsurface flow, the new multilevel estimator achieves gains of more than a factor 60 over standard subset simulation for a practically relevant relative error of 25%. △ Less

Submitted 12 December, 2023; v1 submitted 10 August, 2022; originally announced August 2022.

MSC Class: 65N30; 65C05; 65C40; 35R60

Separating internal and externally-forced contributions to global temperature variability using a Bayesian stochastic energy balance framework

Authors: Maybritt Schillinger, Beatrice Ellerhoff, Robert Scheichl, Kira Rehfeld

Abstract: Earth's temperature variability can be partitioned into internal and externally-forced components. Yet, underlying mechanisms and their relative contributions remain insufficiently understood, especially on decadal to centennial timescales. Important reasons for this are difficulties in isolating internal and externally-forced variability. Here, we provide a physically-motivated emulation of globa… ▽ More Earth's temperature variability can be partitioned into internal and externally-forced components. Yet, underlying mechanisms and their relative contributions remain insufficiently understood, especially on decadal to centennial timescales. Important reasons for this are difficulties in isolating internal and externally-forced variability. Here, we provide a physically-motivated emulation of global mean surface temperature (GMST) variability, which allows for the separation of internal and external variations. To this end, we introduce the ``ClimBayes'' software package, which infers climate parameters from a stochastic energy balance model (EBM) with a Bayesian approach. We apply our method to GMST data from temperature observations and 20 last millennium simulations from climate models of intermediate to high complexity. This yields the best estimates of the EBM's forced and forced + internal response, which we refer to as emulated variability. The timescale-dependent variance is obtained from spectral analysis. In particular, we contrast the emulated forced and forced + internal variance on interannual to centennial timescales with that of the GMST target. Our findings show that a stochastic EBM closely approximates the power spectrum and timescale-dependent variance of GMST as simulated by modern climate models. Small deviations at interannual timescales can be attributed to the simplified representation of internal variability and, in particular, the absence of (pseudo-)oscillatory modes in the stochastic EBM. Altogether, we demonstrate the potential of combining Bayesian inference with conceptual climate models to emulate statistics of climate variables across timescales. △ Less

Submitted 31 January, 2023; v1 submitted 27 June, 2022; originally announced June 2022.

Comments: The following article has been published in Chaos: An Interdisciplinary Journal of Nonlinear Science and can be found at https://aip.scitation.org/doi/abs/10.1063/5.0106123

Journal ref: Chaos 32, 113146 (2022)

Multilevel simulation of hard-sphere mixtures

Authors: Paul B. Rohrbach, Hideki Kobayashi, Robert Scheichl, Nigel B. Wilding, Robert L. Jack

Abstract: We present a multilevel Monte Carlo simulation method for analysing multi-scale physical systems via a hierarchy of coarse-grained representations, to obtain numerically-exact results, at the most detailed level. We apply the method to a mixture of size-asymmetric hard spheres, in the grand canonical ensemble. A three-level version of the method is compared with a previously-studied two-level vers… ▽ More We present a multilevel Monte Carlo simulation method for analysing multi-scale physical systems via a hierarchy of coarse-grained representations, to obtain numerically-exact results, at the most detailed level. We apply the method to a mixture of size-asymmetric hard spheres, in the grand canonical ensemble. A three-level version of the method is compared with a previously-studied two-level version. The extra level interpolates between the full mixture and a coarse-grained description where only the large particles are present -- this is achieved by restricting the small particles to regions close to the large ones. The three-level method improves the performance of the estimator, at fixed computational cost. We analyse the asymptotic variance of the estimator, and discuss the mechanisms for the improved performance. △ Less

Submitted 26 August, 2022; v1 submitted 2 June, 2022; originally announced June 2022.

Comments: 24 pages, authors' accepted version

A Bayesian Approach to Modelling Biological Pattern Formation with Limited Data

Authors: Alexey Kazarnikov, Robert Scheichl, Heikki Haario, Anna Marciniak-Czochra

Abstract: Pattern formation in biological tissues plays an important role in the development of living organisms. Since the classical work of Alan Turing, a pre-eminent way of modelling has been through reaction-diffusion mechanisms. More recently, alternative models have been proposed, that link dynamics of diffusing molecular signals with tissue mechanics. In order to distinguish among different models, t… ▽ More Pattern formation in biological tissues plays an important role in the development of living organisms. Since the classical work of Alan Turing, a pre-eminent way of modelling has been through reaction-diffusion mechanisms. More recently, alternative models have been proposed, that link dynamics of diffusing molecular signals with tissue mechanics. In order to distinguish among different models, they should be compared to experimental observations. However, in many experimental situations only the limiting, stationary regime of the pattern formation process is observable, without knowledge of the transient behaviour or the initial state. The unstable nature of the underlying dynamics in all alternative models seriously complicates model and parameter identification, since small changes in the initial condition lead to distinct stationary patterns. To overcome this problem the initial state of the model can be randomised. In the latter case, fixed values of the model parameters correspond to a family of patterns rather than a fixed stationary solution, and standard approaches to compare pattern data directly with model outputs, e.g., in the least squares sense, are not suitable. Instead, statistical characteristics of the patterns should be compared, which is difficult given the typically limited amount of available data in practical applications. To deal with this problem, we extend a recently developed statistical approach for parameter identification using pattern data, the so-called Correlation Integral Likelihood (CIL) method. We suggest modifications that allow increasing the accuracy of the identification process without resizing the data set. The proposed approach is tested using different classes of pattern formation models. For all considered equations, parallel GPU-based implementations of the numerical solvers with efficient time stepping schemes are provided. △ Less

Submitted 26 January, 2023; v1 submitted 28 March, 2022; originally announced March 2022.

Comments: Minor errors corrected, additional clarification added, results unchanged

Multilevel Delayed Acceptance MCMC

Authors: Mikkel B. Lykkegaard, Tim J. Dodwell, Colin Fox, Grigorios Mingas, Robert Scheichl

Abstract: We develop a novel Markov chain Monte Carlo (MCMC) method that exploits a hierarchy of models of increasing complexity to efficiently generate samples from an unnormalized target distribution. Broadly, the method rewrites the Multilevel MCMC approach of Dodwell et al. (2015) in terms of the Delayed Acceptance (DA) MCMC of Christen & Fox (2005). In particular, DA is extended to use a hierarchy of m… ▽ More We develop a novel Markov chain Monte Carlo (MCMC) method that exploits a hierarchy of models of increasing complexity to efficiently generate samples from an unnormalized target distribution. Broadly, the method rewrites the Multilevel MCMC approach of Dodwell et al. (2015) in terms of the Delayed Acceptance (DA) MCMC of Christen & Fox (2005). In particular, DA is extended to use a hierarchy of models of arbitrary depth, and allow subchains of arbitrary length. We show that the algorithm satisfies detailed balance, hence is ergodic for the target distribution. Furthermore, multilevel variance reduction is derived that exploits the multiple levels and subchains, and an adaptive multilevel correction to coarse-level biases is developed. Three numerical examples of Bayesian inverse problems are presented that demonstrate the advantages of these novel methods. The software and examples are available in PyMC3. △ Less

Submitted 30 August, 2022; v1 submitted 8 February, 2022; originally announced February 2022.

Comments: 29 pages, 12 figures

MSC Class: 62F15; 62M05; 65C05; 65C40

Wavenumber explicit convergence of a multiscale generalized finite element method for heterogeneous Helmholtz problems

Authors: Chupeng Ma, Christian Alber, Robert Scheichl

Abstract: In this paper, a generalized finite element method (GFEM) with optimal local approximation spaces for solving high-frequency heterogeneous Helmholtz problems is systematically studied. The local spaces are built from selected eigenvectors of carefully designed local eigenvalue problems defined on generalized harmonic spaces. At both continuous and discrete levels, $(i)$ wavenumber explicit and nea… ▽ More In this paper, a generalized finite element method (GFEM) with optimal local approximation spaces for solving high-frequency heterogeneous Helmholtz problems is systematically studied. The local spaces are built from selected eigenvectors of carefully designed local eigenvalue problems defined on generalized harmonic spaces. At both continuous and discrete levels, $(i)$ wavenumber explicit and nearly exponential decay rates for local and global approximation errors are obtained without any assumption on the size of subdomains; $(ii)$ a quasi-optimal convergence of the method is established by assuming that the size of subdomains is $O(1/k)$ ($k$ is the wavenumber). A novel resonance effect between the wavenumber and the dimension of local spaces on the decay of error with respect to the oversampling size is implied by the analysis. Furthermore, for fixed dimensions of local spaces, the discrete local errors are proved to converge as $h\rightarrow 0$ ($h$ denoting the mesh size) towards the continuous local errors. The method at the continuous level extends the plane wave partition of unity method [I. Babuska and J. M. Melenk, Int.\;J.\;Numer.\;Methods Eng., 40 (1997), pp.~727--758] to the heterogeneous-coefficients case, and at the discrete level, it delivers an efficient non-iterative domain decomposition method for solving discrete Helmholtz problems resulting from standard FE discretizations. Numerical results are provided to confirm the theoretical analysis and to validate the proposed method. △ Less

Submitted 14 September, 2022; v1 submitted 20 December, 2021; originally announced December 2021.

Overlapping Schwarz methods with GenEO coarse spaces for indefinite and non-self-adjoint problems

Authors: Niall Bootland, Victorita Dolean, Ivan G. Graham, Chupeng Ma, Robert Scheichl

Abstract: GenEO (`Generalised Eigenvalue problems on the Overlap') is a method for computing an operator-dependent spectral coarse space to be combined with local solves on subdomains to form a robust parallel domain decomposition preconditioner for elliptic PDEs. It has previously been proved, in the self-adjoint and positive-definite case, that this method, when used as a preconditioner for conjugate grad… ▽ More GenEO (`Generalised Eigenvalue problems on the Overlap') is a method for computing an operator-dependent spectral coarse space to be combined with local solves on subdomains to form a robust parallel domain decomposition preconditioner for elliptic PDEs. It has previously been proved, in the self-adjoint and positive-definite case, that this method, when used as a preconditioner for conjugate gradients, yields iteration numbers which are completely independent of the heterogeneity of the coefficient field of the partial differential operator. We extend this theory to the case of convection-diffusion-reaction problems, which may be non-self-adjoint and indefinite, and whose discretisations are solved with preconditioned GMRES. The GenEO coarse space is defined here using a generalised eigenvalue problem based on a self-adjoint and positive-definite subproblem. We prove estimates on GMRES iteration counts which are independent of the variation of the coefficient of the diffusion term in the operator and depend only very mildly on variations of the other coefficients. These are proved under the assumption that the subdomain diameter is sufficiently small and the eigenvalue tolerance for building the coarse space is sufficiently large. While the iteration number estimates do grow as the non-self-adjointness and indefiniteness of the operator increases, practical tests indicate the deterioration is much milder. Thus we obtain an iterative solver which is efficient in parallel and very effective for a wide range of convection--diffusion--reaction problems. △ Less

Submitted 25 March, 2022; v1 submitted 26 October, 2021; originally announced October 2021.

MSC Class: 65N22; 65N55; 65F10

Journal ref: IMA J. Numer. Anal., 43(4), 1899-1936 (2023)

High Performance Uncertainty Quantification with Parallelized Multilevel Markov Chain Monte Carlo

Authors: Linus Seelinger, Anne Reinarz, Leonhard Rannabauer, Michael Bader, Peter Bastian, Robert Scheichl

Abstract: Numerical models of complex real-world phenomena often necessitate High Performance Computing (HPC). Uncertainties increase problem dimensionality further and pose even greater challenges. We present a parallelization strategy for multilevel Markov chain Monte Carlo, a state-of-the-art, algorithmically scalable Uncertainty Quantification (UQ) algorithm for Bayesian inverse problems, and a new so… ▽ More Numerical models of complex real-world phenomena often necessitate High Performance Computing (HPC). Uncertainties increase problem dimensionality further and pose even greater challenges. We present a parallelization strategy for multilevel Markov chain Monte Carlo, a state-of-the-art, algorithmically scalable Uncertainty Quantification (UQ) algorithm for Bayesian inverse problems, and a new software framework allowing for large-scale parallelism across forward model evaluations and the UQ algorithms themselves. The main scalability challenge presents itself in the form of strong data dependencies introduced by the MLMCMC method, prohibiting trivial parallelization. Our software is released as part of the modular and open-source MIT UQ Library (MUQ), and can easily be coupled with arbitrary user codes. We demonstrate it using the DUNE and the ExaHyPE Engine. The latter provides a realistic, large-scale tsunami model in which identify the source of a tsunami from buoy-elevation data. △ Less

Submitted 30 July, 2021; originally announced July 2021.

Error estimates for fully discrete generalized FEMs with locally optimal spectral approximations

Authors: Chupeng Ma, Robert Scheichl

Abstract: This paper is concerned with error estimates of the fully discrete generalized finite element method (GFEM) with optimal local approximation spaces for solving elliptic problems with heterogeneous coefficients. The local approximation spaces are constructed using eigenvectors of local eigenvalue problems solved by the finite element method on some sufficiently fine mesh with mesh size $h$. The err… ▽ More This paper is concerned with error estimates of the fully discrete generalized finite element method (GFEM) with optimal local approximation spaces for solving elliptic problems with heterogeneous coefficients. The local approximation spaces are constructed using eigenvectors of local eigenvalue problems solved by the finite element method on some sufficiently fine mesh with mesh size $h$. The error bound of the discrete GFEM approximation is proved to converge as $h\rightarrow 0$ towards that of the continuous GFEM approximation, which was shown to decay nearly exponentially in previous works. Moreover, even for fixed mesh size $h$, a nearly exponential rate of convergence of the local approximation errors with respect to the dimension of the local spaces is established. An efficient and accurate method for solving the discrete eigenvalue problems is proposed by incorporating the discrete $A$-harmonic constraint directly into the eigensolver. Numerical experiments are carried out to confirm the theoretical results and to demonstrate the effectiveness of the method. △ Less

Submitted 30 September, 2021; v1 submitted 21 July, 2021; originally announced July 2021.

MSC Class: 65M60; 65N15; 65N55

Critical point for de-mixing of binary hard spheres

Authors: Hideki Kobayashi, Paul B. Rohrbach, Robert Scheichl, Nigel B. Wilding, Robert L. Jack

Abstract: We use a two-level simulation method to analyse the critical point associated with demixing of binary hard sphere mixtures. The method exploits an accurate coarse-grained model with two-body and three-body effective interactions. Using this model within the two-level methodology allows computation of properties of the full (fine-grained) mixture. The critical point is located by computing the prob… ▽ More We use a two-level simulation method to analyse the critical point associated with demixing of binary hard sphere mixtures. The method exploits an accurate coarse-grained model with two-body and three-body effective interactions. Using this model within the two-level methodology allows computation of properties of the full (fine-grained) mixture. The critical point is located by computing the probability distribution for the number of large particles in the grand canonical ensemble, and matching to the universal form for the $3d$ Ising universality class. The results have a strong and unexpected dependence on the size ratio between large and small particles, which is related to three-body effective interactions, and the geometry of the underlying hard sphere packings. △ Less

Submitted 2 July, 2021; originally announced July 2021.

Comments: 11 pages

Journal ref: Phys. Rev. E 104, 044603 (2021)

Multilevel Spectral Domain Decomposition

Authors: Peter Bastian, Robert Scheichl, Linus Seelinger, Arne Strehlow

Abstract: Highly heterogeneous, anisotropic coefficients, e.g. in the simulation of carbon-fibre composite components, can lead to extremely challenging finite element systems. Direct solvers for the resulting large and sparse linear systems suffer from severe memory requirements and limited parallel scalability, while iterative solvers in general lack robustness. Two-level spectral domain decomposition met… ▽ More Highly heterogeneous, anisotropic coefficients, e.g. in the simulation of carbon-fibre composite components, can lead to extremely challenging finite element systems. Direct solvers for the resulting large and sparse linear systems suffer from severe memory requirements and limited parallel scalability, while iterative solvers in general lack robustness. Two-level spectral domain decomposition methods can provide such robustness for symmetric positive definite linear systems, by using coarse spaces based on independent generalized eigenproblems in the subdomains. Rigorous condition number bounds are independent of mesh size, number of subdomains, as well as coefficient contrast. However, their parallel scalability is still limited by the fact that (in order to guarantee robustness) the coarse problem is solved via a direct method. In this paper, we introduce a multilevel variant in the context of subspace correction methods and provide a general convergence theory for its robust convergence for abstract, elliptic variational problems. Assumptions of the theory are verified for conforming, as well as for discontinuous Galerkin methods applied to a scalar diffusion problem. Numerical results illustrate the performance of the method for two- and three-dimensional problems and for various discretization schemes, in the context of scalar diffusion and linear elasticity. △ Less

Submitted 15 June, 2021; v1 submitted 11 June, 2021; originally announced June 2021.

Comments: 25 pages

MSC Class: 65N55 (Primary) 65F08; 65F10 (Secondary)

GenEO coarse spaces for heterogeneous indefinite elliptic problems

Authors: Niall Bootland, Victorita Dolean, Ivan G. Graham, Chupeng Ma, Robert Scheichl

Abstract: Motivated by recent work on coarse spaces for Helmholtz problems, we provide in this paper a comparative study on the use of spectral coarse spaces of GenEO type for heterogeneous indefinite elliptic problems within an additive overlapping Schwarz method. In particular, we focus here on two different but related formulations of local generalised eigenvalue problems and compare their performance nu… ▽ More Motivated by recent work on coarse spaces for Helmholtz problems, we provide in this paper a comparative study on the use of spectral coarse spaces of GenEO type for heterogeneous indefinite elliptic problems within an additive overlapping Schwarz method. In particular, we focus here on two different but related formulations of local generalised eigenvalue problems and compare their performance numerically. Even though their behaviour seems to be very similar for several well-known heterogeneous test cases that are mildly indefinite, only one of the coarse spaces has so far been analysed theoretically, while the other one leads to a significantly more robust domain decomposition method when the indefiniteness is increased. We present a summary of upcoming results developing such a theory and describe how the numerical experiments illustrate it. △ Less

Submitted 7 July, 2021; v1 submitted 30 March, 2021; originally announced March 2021.

Comments: 26th International Domain Decomposition Conference (DD26) proceedings

MSC Class: 65N55; 65N25; 65F08; 65F10; 65-06

Journal ref: Brenner, S.C., Chung, E., Klawonn, A., Kwok, F., Xu, J., Zou, J. (Eds.) Domain Decomposition Methods in Science and Engineering XXVI. Lecture Notes in Computational Science and Engineering, Vol. 145. Springer, Cham, pp. 117-125 (2022)

Novel design and analysis of generalized FE methods based on locally optimal spectral approximations

Authors: Chupeng Ma, Robert Scheichl, Tim Dodwell

Abstract: In this paper, the generalized finite element method (GFEM) for solving second order elliptic equations with rough coefficients is studied. New optimal local approximation spaces for GFEMs based on local eigenvalue problems involving a partition of unity are presented. These new spaces have advantages over those proposed in [I. Babuska and R. Lipton, Multiscale Model.\;\,Simul., 9 (2011), pp.~373-… ▽ More In this paper, the generalized finite element method (GFEM) for solving second order elliptic equations with rough coefficients is studied. New optimal local approximation spaces for GFEMs based on local eigenvalue problems involving a partition of unity are presented. These new spaces have advantages over those proposed in [I. Babuska and R. Lipton, Multiscale Model.\;\,Simul., 9 (2011), pp.~373--406]. First, in addition to a nearly exponential decay rate of the local approximation errors with respect to the dimensions of the local spaces, the rate of convergence with respect to the size of the oversampling region is also established. Second, the theoretical results hold for problems with mixed boundary conditions defined on general Lipschitz domains. Finally, an efficient and easy-to-implement technique for generating the discrete $A$-harmonic spaces is proposed which relies on solving an eigenvalue problem associated with the Dirichlet-to-Neumann operator, leading to a substantial reduction in computational cost. Numerical experiments are presented to support the theoretical analysis and to confirm the effectiveness of the new method. △ Less

Submitted 21 December, 2021; v1 submitted 17 March, 2021; originally announced March 2021.

Multilevel quasi-Monte Carlo for random elliptic eigenvalue problems II: Efficient algorithms and numerical results

Authors: Alexander D. Gilbert, Robert Scheichl

Abstract: Stochastic PDE eigenvalue problems often arise in the field of uncertainty quantification, whereby one seeks to quantify the uncertainty in an eigenvalue, or its eigenfunction. In this paper we present an efficient multilevel quasi-Monte Carlo (MLQMC) algorithm for computing the expectation of the smallest eigenvalue of an elliptic eigenvalue problem with stochastic coefficients. Each sample evalu… ▽ More Stochastic PDE eigenvalue problems often arise in the field of uncertainty quantification, whereby one seeks to quantify the uncertainty in an eigenvalue, or its eigenfunction. In this paper we present an efficient multilevel quasi-Monte Carlo (MLQMC) algorithm for computing the expectation of the smallest eigenvalue of an elliptic eigenvalue problem with stochastic coefficients. Each sample evaluation requires the solution of a PDE eigenvalue problem, and so tackling this problem in practice is notoriously computationally difficult. We speed up the approximation of this expectation in four ways: we use a multilevel variance reduction scheme to spread the work over a hierarchy of FE meshes and truncation dimensions; we use QMC methods to efficiently compute the expectations on each level; we exploit the smoothness in parameter space and reuse the eigenvector from a nearby QMC point to reduce the number of iterations of the eigensolver; and we utilise a two-grid discretisation scheme to obtain the eigenvalue on the fine mesh with a single linear solve. The full error analysis of a basic MLQMC algorithm is given in the companion paper [Gilbert and Scheichl, 2022], and so in this paper we focus on how to further improve the efficiency and provide theoretical justification for using nearby QMC points and two-grid methods. Numerical results are presented that show the efficiency of our algorithm, and also show that the four strategies we employ are complementary. △ Less

Submitted 6 October, 2022; v1 submitted 4 March, 2021; originally announced March 2021.

Multilevel Delayed Acceptance MCMC with an Adaptive Error Model in PyMC3

Authors: Mikkel B. Lykkegaard, Grigorios Mingas, Robert Scheichl, Colin Fox, Tim J. Dodwell

Abstract: Uncertainty Quantification through Markov Chain Monte Carlo (MCMC) can be prohibitively expensive for target probability densities with expensive likelihood functions, for instance when the evaluation it involves solving a Partial Differential Equation (PDE), as is the case in a wide range of engineering applications. Multilevel Delayed Acceptance (MLDA) with an Adaptive Error Model (AEM) is a nov… ▽ More Uncertainty Quantification through Markov Chain Monte Carlo (MCMC) can be prohibitively expensive for target probability densities with expensive likelihood functions, for instance when the evaluation it involves solving a Partial Differential Equation (PDE), as is the case in a wide range of engineering applications. Multilevel Delayed Acceptance (MLDA) with an Adaptive Error Model (AEM) is a novel approach, which alleviates this problem by exploiting a hierarchy of models, with increasing complexity and cost, and correcting the inexpensive models on-the-fly. The method has been integrated within the open-source probabilistic programming package PyMC3 and is available in the latest development version. In this paper, the algorithm is presented along with an illustrative example. △ Less

Submitted 10 December, 2020; originally announced December 2020.

Comments: 8 pages, 4 figures, accepted for Machine Learning for Engineering Modeling, Simulation, and Design Workshop at Neural Information Processing Systems 2020

Multilevel quasi-Monte Carlo for random elliptic eigenvalue problems I: Regularity and error analysis

Authors: Alexander D. Gilbert, Robert Scheichl

Abstract: Stochastic PDE eigenvalue problems are useful models for quantifying the uncertainty in several applications from the physical sciences and engineering, e.g., structural vibration analysis, the criticality of a nuclear reactor or photonic crystal structures. In this paper we present a multilevel quasi-Monte Carlo (MLQMC) method for approximating the expectation of the minimal eigenvalue of an elli… ▽ More Stochastic PDE eigenvalue problems are useful models for quantifying the uncertainty in several applications from the physical sciences and engineering, e.g., structural vibration analysis, the criticality of a nuclear reactor or photonic crystal structures. In this paper we present a multilevel quasi-Monte Carlo (MLQMC) method for approximating the expectation of the minimal eigenvalue of an elliptic eigenvalue problem with coefficients that are given as a series expansion of countably-many stochastic parameters. The MLQMC algorithm is based on a hierarchy of discretisations of the spatial domain and truncations of the dimension of the stochastic parameter domain. To approximate the expectations, randomly shifted lattice rules are employed. This paper is primarily dedicated to giving a rigorous analysis of the error of this algorithm. A key step in the error analysis requires bounds on the mixed derivatives of the eigenfunction with respect to both the stochastic and spatial variables simultaneously. Under stronger smoothness assumptions on the parametric dependence, our analysis also extends to multilevel higher-order quasi-Monte Carlo rules. An accompanying paper [Gilbert and Scheichl, 2022], focusses on practical extensions of the MLQMC algorithm to improve efficiency, and presents numerical results. △ Less

Submitted 6 October, 2022; v1 submitted 2 October, 2020; originally announced October 2020.

Multilevel Monte Carlo for quantum mechanics on a lattice

Authors: Karl Jansen, Eike Hermann Müller, Robert Scheichl

Abstract: Monte Carlo simulations of quantum field theories on a lattice become increasingly expensive as the continuum limit is approached since the cost per independent sample grows with a high power of the inverse lattice spacing. Simulations on fine lattices suffer from critical slowdown, the rapid growth of autocorrelations in the Markov chain. This causes a strong increase in the number of lattice con… ▽ More Monte Carlo simulations of quantum field theories on a lattice become increasingly expensive as the continuum limit is approached since the cost per independent sample grows with a high power of the inverse lattice spacing. Simulations on fine lattices suffer from critical slowdown, the rapid growth of autocorrelations in the Markov chain. This causes a strong increase in the number of lattice configurations that have to be generated to obtain statistically significant results. This paper discusses hierarchical sampling methods to tame the growth in autocorrelations. Combined with multilevel variance reduction, this significantly reduces the computational cost of simulations for given tolerances $εいぷしろん_{\text{disc}}$ on the discretisation error and $εいぷしろん_{\text{stat}}$ on the statistical error. For observables with lattice errors of order $αあるふぁ$ and integrated autocorrelation times that grow like $τたう_{\mathrm{int}}\propto a^{-z}$, multilevel Monte Carlo (MLMC) reduces the cost from $\mathcal{O}(εいぷしろん_{\text{stat}}^{-2}εいぷしろん_{\text{disc}}^{-(1+z)/αあるふぁ})$ to $\mathcal{O}(εいぷしろん_{\text{stat}}^{-2}\vert\log εいぷしろん_{\text{disc}} \vert^2+εいぷしろん_{\text{disc}}^{-1/αあるふぁ})$ or $\mathcal{O}(εいぷしろん_{\text{stat}}^{-2}+εいぷしろん_{\text{disc}}^{-1/αあるふぁ})$. Higher gains are expected for simulations of quantum field theories in $D$ dimensions. The efficiency of the approach is demonstrated on two model systems, including a topological oscillator that is badly affected by critical slowdown from topological charge freezing. On fine lattices, the new methods are orders of magnitude faster than standard Hybrid Monte Carlo sampling. For high resolutions, MLMC can be used to accelerate even the cluster algorithm for the topological oscillator. Performance is further improved through perturbative matching which guarantees efficient coupling of theories on the multilevel hierarchy. △ Less

Submitted 17 November, 2020; v1 submitted 7 August, 2020; originally announced August 2020.

Comments: 23 pages, 17 figures, 5 tables

MSC Class: 81-08; 81T25; 65Y20; 60J22 ACM Class: F.2; J.2

Journal ref: Phys. Rev. D 102, 114512 (2020)

Rank Bounds for Approximating Gaussian Densities in the Tensor-Train Format

Authors: Paul B. Rohrbach, Sergey Dolgov, Lars Grasedyck, Robert Scheichl

Abstract: Low-rank tensor approximations have shown great potential for uncertainty quantification in high dimensions, for example, to build surrogate models that can be used to speed up large-scale inference problems (Eigel et al., Inverse Problems 34, 2018; Dolgov et al., Statistics & Computing 30, 2020). The feasibility and efficiency of such approaches depends critically on the rank that is necessary to… ▽ More Low-rank tensor approximations have shown great potential for uncertainty quantification in high dimensions, for example, to build surrogate models that can be used to speed up large-scale inference problems (Eigel et al., Inverse Problems 34, 2018; Dolgov et al., Statistics & Computing 30, 2020). The feasibility and efficiency of such approaches depends critically on the rank that is necessary to represent or approximate the underlying distribution. In this paper, a-priori rank bounds for approximations in the functional tensor-train representation for the case of Gaussian models are developed. It is shown that under suitable conditions on the precision matrix, the Gaussian density can be approximated to high accuracy without suffering from an exponential growth of complexity as the dimension increases. These results provide a rigorous justification of the suitability and the limitations of low-rank tensor methods in a simple but important model case. Numerical experiments confirm that the rank bounds capture the qualitative behavior of the rank structure when varying the parameters of the precision matrix and the accuracy of the approximation. Finally, the practical relevance of the theoretical results is demonstrated in the context of a Bayesian filtering problem. △ Less

Submitted 27 November, 2020; v1 submitted 22 January, 2020; originally announced January 2020.

Comments: 25 pages, 6 figures

MSC Class: 15A23; 15A69; 65C60; 65D32; 65D15; 41A10

Multilevel Dimension-Independent Likelihood-Informed MCMC for Large-Scale Inverse Problems

Authors: Tiangang Cui, Gianluca Detommaso, Robert Scheichl

Abstract: We present a non-trivial integration of dimension-independent likelihood-informed (DILI) MCMC (Cui, Law, Marzouk, 2016) and the multilevel MCMC (Dodwell et al., 2015) to explore the hierarchy of posterior distributions. This integration offers several advantages: First, DILI-MCMC employs an intrinsic likelihood-informed subspace (LIS) (Cui et al., 2014) -- which involves a number of forward and ad… ▽ More We present a non-trivial integration of dimension-independent likelihood-informed (DILI) MCMC (Cui, Law, Marzouk, 2016) and the multilevel MCMC (Dodwell et al., 2015) to explore the hierarchy of posterior distributions. This integration offers several advantages: First, DILI-MCMC employs an intrinsic likelihood-informed subspace (LIS) (Cui et al., 2014) -- which involves a number of forward and adjoint model simulations -- to design accelerated operator-weighted proposals. By exploiting the multilevel structure of the discretised parameters and discretised forward models, we design a Rayleigh-Ritz procedure to significantly reduce the computational effort in building the LIS and operating with DILI proposals. Second, the resulting DILI-MCMC can drastically improve the sampling efficiency of MCMC at each level, and hence reduce the integration error of the multilevel algorithm for fixed CPU time. Numerical results confirm the improved computational efficiency of the multilevel DILI approach. △ Less

Submitted 29 November, 2023; v1 submitted 28 October, 2019; originally announced October 2019.

Multilevel Monte Carlo Simulations of Composite Structures with Uncertain Manufacturing Defects

Authors: T. J. Dodwell, S. Kinston, R. Butler, R. T. Haftka, Nam H. Kim, R. Scheichl

Abstract: By adopting a Multilevel Monte Carlo (MLMC) framework, we show that only a handful of costly fine scale computations are needed to accurately estimate statistics of the failure of a composite structure, as opposed to the thousands typically needed in classical Monte Carlo analyses. We introduce the MLMC method, compare its theoretical complexity with classical Monte Carlo, and give a simple-to-imp… ▽ More By adopting a Multilevel Monte Carlo (MLMC) framework, we show that only a handful of costly fine scale computations are needed to accurately estimate statistics of the failure of a composite structure, as opposed to the thousands typically needed in classical Monte Carlo analyses. We introduce the MLMC method, compare its theoretical complexity with classical Monte Carlo, and give a simple-to-implement algorithm which includes a simple extension called MLMC with selective refinement to efficiently calculated structural failure probabilities. To demonstrate the huge computational gains we present two benchmark problems in composites: (1) the effects of fibre waviness on the compressive strength of a composite material, (2) uncertain buckling performance of a composite panel with uncertain ply orientations. For our most challenging test case, estimating a rare ($\sim 1/150$) probability of buckling failure of a composite panel, we see a speed-up factor $> 1000$. Our approach distributed over $1024$ processors reduces the computation time from $218$ days to just $4.5$ hours. This level of speed up makes stochastic simulations that would otherwise be unthinkable now possible.? △ Less

Submitted 24 July, 2019; originally announced July 2019.

Comments: 23 Pages, 7 Figures

Correction of coarse-graining errors by a two-level method: application to the Asakura-Oosawa model

Authors: Hideki Kobayashi, Paul B. Rohrbach, Robert Scheichl, Nigel B. Wilding, Robert L. Jack

Abstract: We present a method that exploits self-consistent simulation of coarse-grained and fine-grained models, in order to analyse properties of physical systems. The method uses the coarse-grained model to obtain a first estimate of the quantity of interest, before computing a correction by analysing properties of the fine system. We illustrate the method by applying it to the Asakura-Oosawa (AO) model… ▽ More We present a method that exploits self-consistent simulation of coarse-grained and fine-grained models, in order to analyse properties of physical systems. The method uses the coarse-grained model to obtain a first estimate of the quantity of interest, before computing a correction by analysing properties of the fine system. We illustrate the method by applying it to the Asakura-Oosawa (AO) model of colloid-polymer mixtures. We show that the liquid-vapour critical point in that system is affected by three-body interactions which are neglected in the corresponding coarse-grained model. We analyse the size of this effect and the nature of the three-body interactions. We also analyse the accuracy of the method, as a function of the associated computational effort. △ Less

Submitted 18 July, 2019; originally announced July 2019.

Comments: 17 pages

A High-Performance Implementation of a Robust Preconditioner for Heterogeneous Problems

Authors: Linus Seelinger, Anne Reinarz, Robert Scheichl

Abstract: We present an efficient implementation of the highly robust and scalable GenEO preconditioner in the high-performance PDE framework DUNE. The GenEO coarse space is constructed by combining low energy solutions of a local generalised eigenproblem using a partition of unity. In this paper we demonstrate both weak and strong scaling for the GenEO solver on over 15,000 cores by solving an industrially… ▽ More We present an efficient implementation of the highly robust and scalable GenEO preconditioner in the high-performance PDE framework DUNE. The GenEO coarse space is constructed by combining low energy solutions of a local generalised eigenproblem using a partition of unity. In this paper we demonstrate both weak and strong scaling for the GenEO solver on over 15,000 cores by solving an industrially motivated problem with over 200 million degrees of freedom. Further, we show that for highly complex parameter distributions arising in certain real-world applications, established methods become intractable while GenEO remains fully effective. The purpose of this paper is two-fold: to demonstrate the robustness and high parallel efficiency of the solver and to document the technical details that are crucial to the efficiency of the code. △ Less

Submitted 16 June, 2020; v1 submitted 26 June, 2019; originally announced June 2019.

Comments: Fixed error in definition of partition of unity

MSC Class: 68U20

Unified Analysis of Periodization-Based Sampling Methods for Matérn Covariances

Authors: Markus Bachmayr, Ivan G. Graham, Van Kien Nguyen, Robert Scheichl

Abstract: The periodization of a stationary Gaussian random field on a sufficiently large torus comprising the spatial domain of interest is the basis of various efficient computational methods, such as the classical circulant embedding technique using the fast Fourier transform for generating samples on uniform grids. For the family of Matérn covariances with smoothness index $νにゅー$ and correlation length… ▽ More The periodization of a stationary Gaussian random field on a sufficiently large torus comprising the spatial domain of interest is the basis of various efficient computational methods, such as the classical circulant embedding technique using the fast Fourier transform for generating samples on uniform grids. For the family of Matérn covariances with smoothness index $νにゅー$ and correlation length $λらむだ$, we analyse the nonsmooth periodization (corresponding to classical circulant embedding) and an alternative procedure using a smooth truncation of the covariance function. We solve two open problems: the first concerning the $νにゅー$-dependent asymptotic decay of eigenvalues of the resulting circulant in the nonsmooth case, the second concerning the required size in terms of $νにゅー$, $λらむだ$ of the torus when using a smooth periodization. In doing this we arrive at a complete characterisation of the performance of these two approaches. Both our theoretical estimates and the numerical tests provided here show substantial advantages of smooth truncation. △ Less

Submitted 25 August, 2020; v1 submitted 31 May, 2019; originally announced May 2019.

Comments: 24 pages, 3 figures

MSC Class: 60G15; 60G60; 42B05; 65T40

HINT: Hierarchical Invertible Neural Transport for Density Estimation and Bayesian Inference

Authors: Jakob Kruse, Gianluca Detommaso, Ullrich Köthe, Robert Scheichl

Abstract: Many recent invertible neural architectures are based on coupling block designs where variables are divided in two subsets which serve as inputs of an easily invertible (usually affine) triangular transformation. While such a transformation is invertible, its Jacobian is very sparse and thus may lack expressiveness. This work presents a simple remedy by noting that subdivision and (affine) couplin… ▽ More Many recent invertible neural architectures are based on coupling block designs where variables are divided in two subsets which serve as inputs of an easily invertible (usually affine) triangular transformation. While such a transformation is invertible, its Jacobian is very sparse and thus may lack expressiveness. This work presents a simple remedy by noting that subdivision and (affine) coupling can be repeated recursively within the resulting subsets, leading to an efficiently invertible block with dense, triangular Jacobian. By formulating our recursive coupling scheme via a hierarchical architecture, HINT allows sampling from a joint distribution p(y,x) and the corresponding posterior p(x|y) using a single invertible network. We evaluate our method on some standard data sets and benchmark its full power for density estimation and Bayesian inference on a novel data set of 2D shapes in Fourier parameterization, which enables consistent visualization of samples for different dimensionalities. △ Less

Submitted 25 May, 2021; v1 submitted 25 May, 2019; originally announced May 2019.

Comments: Published at AAAI 2021

Full Error Analysis and Uncertainty Quantification for the Heterogeneous Transport Equation in Slab Geometry

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

Abstract: We present an analysis of multilevel Monte Carlo techniques for the forward problem of uncertainty quantification for the radiative transport equation, when the coefficients ({\em cross-sections}) are heterogenous random fields. To do this, we first give a new error analysis for the combined spatial and angular discretisation in the deterministic case, with error estimates which are explicit in th… ▽ More We present an analysis of multilevel Monte Carlo techniques for the forward problem of uncertainty quantification for the radiative transport equation, when the coefficients ({\em cross-sections}) are heterogenous random fields. To do this, we first give a new error analysis for the combined spatial and angular discretisation in the deterministic case, with error estimates which are explicit in the coefficients (and allow for very low regularity and jumps). This detailed error analysis is done for the 1D space - 1D angle slab geometry case with classical diamond differencing. Under reasonable assumptions on the statistics of the coefficients, we then prove an error estimate for the random problem in a suitable Bochner space. Because the problem is not self-adjoint, stability can only be proved under a path-dependent mesh resolution condition. This means that, while the Bochner space error estimate is of order $\mathcal{O}(h^ηいーた)$ for some $ηいーた$, where $h$ is a (deterministically chosen) mesh diameter, smaller mesh sizes might be needed for some realisations. Under reasonable assumptions we show that the expected cost for computing a typical quantity of interest remains of the same order as for a single sample. This leads to rigorous complexity estimates for Monte Carlo and multilevel Monte Carlo: For particular linear solvers, the multilevel version gives up to two orders of magnitude improvement over Monte Carlo. We provide numerical results supporting the theory. △ Less

Submitted 14 January, 2020; v1 submitted 28 March, 2019; originally announced March 2019.

MSC Class: 65N12; 65R99; 65C30; 65C05

A fully adaptive multilevel stochastic collocation strategy for solving elliptic PDEs with random data

Authors: Jens Lang, Robert Scheichl, David Silvester

Abstract: We propose and analyse a fully adaptive strategy for solving elliptic PDEs with random data in this work. A hierarchical sequence of adaptive mesh refinements for the spatial approximation is combined with adaptive anisotropic sparse Smolyak grids in the stochastic space in such a way as to minimize the computational cost. The novel aspect of our strategy is that the hierarchy of spatial approxima… ▽ More We propose and analyse a fully adaptive strategy for solving elliptic PDEs with random data in this work. A hierarchical sequence of adaptive mesh refinements for the spatial approximation is combined with adaptive anisotropic sparse Smolyak grids in the stochastic space in such a way as to minimize the computational cost. The novel aspect of our strategy is that the hierarchy of spatial approximations is sample dependent so that the computational effort at each collocation point can be optimised individually. We outline a rigorous analysis for the convergence and computational complexity of the adaptive multilevel algorithm and we provide optimal choices for error tolerances at each level. Two numerical examples demonstrate the reliability of the error control and the significant decrease in the complexity that arises when compared to single level algorithms and multilevel algorithms that employ adaptivity solely in the spatial discretisation or in the collocation procedure. △ Less

Submitted 4 December, 2019; v1 submitted 9 February, 2019; originally announced February 2019.

Comments: 26 pages, 7 figures

MSC Class: 65C20; 65C30; 65N35; 65M75

Bounding the spectral gap for an elliptic eigenvalue problem with uniformly bounded stochastic coefficients

Authors: Alexander D. Gilbert, Ivan G. Graham, Robert Scheichl, Ian H. Sloan

Abstract: A key quantity that occurs in the error analysis of several numerical methods for eigenvalue problems is the distance between the eigenvalue of interest and the next nearest eigenvalue. When we are interested in the smallest or fundamental eigenvalue, we call this the spectral or fundamental gap. In a recent manuscript [Gilbert et al., arXiv:1808.02639], the current authors, together with Frances… ▽ More A key quantity that occurs in the error analysis of several numerical methods for eigenvalue problems is the distance between the eigenvalue of interest and the next nearest eigenvalue. When we are interested in the smallest or fundamental eigenvalue, we call this the spectral or fundamental gap. In a recent manuscript [Gilbert et al., arXiv:1808.02639], the current authors, together with Frances Kuo, studied an elliptic eigenvalue problem with homogeneous Dirichlet boundary conditions, and with coefficients that depend on an infinite number of uniformly distributed stochastic parameters. In this setting, the eigenvalues, and in turn the eigenvalue gap, also depend on the stochastic parameters. Hence, for a robust error analysis one needs to be able to bound the gap over all possible realisations of the parameters, and because the gap depends on infinitely-many random parameters, this is not trivial. This short note presents, in a simplified setting, an important result that was shown in the paper above. Namely, that, under certain decay assumptions on the coefficient, the spectral gap of such a random elliptic eigenvalue problem can be bounded away from 0, uniformly over the entire infinite-dimensional parameter space. △ Less

Submitted 29 January, 2019; originally announced January 2019.

Comments: 12 pages, 5 figures. arXiv admin note: text overlap with arXiv:1808.02639

High-performance dune modules for solving large-scale, strongly anisotropic elliptic problems with applications to aerospace composites

Authors: Richard Butler, Tim Dodwell, Anne Reinarz, Anhad Sandhu, Robert Scheichl, Linus Seelinger

Abstract: The key innovation in this paper is an open-source, high-performance iterative solver for high contrast, strongly anisotropic elliptic partial differential equations implemented within dune-pdelab. The iterative solver exploits a robust, scalable two-level additive Schwarz preconditioner, GenEO (Spillane et al. 2014). The development of this solver has been motivated by the need to overcome the li… ▽ More The key innovation in this paper is an open-source, high-performance iterative solver for high contrast, strongly anisotropic elliptic partial differential equations implemented within dune-pdelab. The iterative solver exploits a robust, scalable two-level additive Schwarz preconditioner, GenEO (Spillane et al. 2014). The development of this solver has been motivated by the need to overcome the limitations of commercially available modeling tools for solving structural analysis simulations in aerospace composite applications. Our software toolbox dune-composites encapsulates the mathematical complexities of the underlying packages within an efficient C++ framework, providing an application interface to our new high-performance solver. We illustrate its use on a range of industrially motivated examples, which should enable other scientists to build on and extend dune-composites and the GenEO preconditioner for use in their own applications. We demonstrate the scalability of the solver on more than 15,000 cores of the UK national supercomputer Archer, solving an aerospace composite problem with over 200 million degrees of freedom in a few minutes. This scale of computation brings composites problems that would otherwise be unthinkable into the feasible range. To demonstrate the wider applicability of the new solver, we also confirm the robustness and scalability of the solver on SPE10, a challenging benchmark in subsurface flow/reservoir simulation. △ Less

Submitted 29 October, 2019; v1 submitted 16 January, 2019; originally announced January 2019.

Approximation and sampling of multivariate probability distributions in the tensor train decomposition

Authors: Sergey Dolgov, Karim Anaya-Izquierdo, Colin Fox, Robert Scheichl

Abstract: General multivariate distributions are notoriously expensive to sample from, particularly the high-dimensional posterior distributions in PDE-constrained inverse problems. This paper develops a sampler for arbitrary continuous multivariate distributions that is based on low-rank surrogates in the tensor-train format. We construct a tensor-train approximation to the target probability density funct… ▽ More General multivariate distributions are notoriously expensive to sample from, particularly the high-dimensional posterior distributions in PDE-constrained inverse problems. This paper develops a sampler for arbitrary continuous multivariate distributions that is based on low-rank surrogates in the tensor-train format. We construct a tensor-train approximation to the target probability density function using the cross interpolation, which requires a small number of function evaluations. For sufficiently smooth distributions the storage required for the TT approximation is moderate, scaling linearly with dimension. The structure of the tensor-train surrogate allows efficient sampling by the conditional distribution method. Unbiased estimates may be calculated by correcting the transformed random seeds using a Metropolis--Hastings accept/reject step. Moreover, one can use a more efficient quasi-Monte Carlo quadrature that may be corrected either by a control-variate strategy, or by importance weighting. We show that the error in the tensor-train approximation propagates linearly into the Metropolis--Hastings rejection rate and the integrated autocorrelation time of the resulting Markov chain. These methods are demonstrated in three computed examples: fitting failure time of shock absorbers; a PDE-constrained inverse diffusion problem; and sampling from the Rosenbrock distribution. The delayed rejection adaptive Metropolis (DRAM) algorithm is used as a benchmark. We find that the importance-weight corrected quasi-Monte Carlo quadrature performs best in all computed examples, and is orders-of-magnitude more efficient than DRAM across a wide range of approximation accuracies and sample sizes. Indeed, all the methods developed here significantly outperform DRAM in all computed examples. △ Less

Submitted 3 July, 2019; v1 submitted 2 October, 2018; originally announced October 2018.

Comments: 32 pages

MSC Class: 65D15; 65D32; 65C05; 65C40; 65C60; 62F15; 15A69; 15A23

Analysis of quasi-Monte Carlo methods for elliptic eigenvalue problems with stochastic coefficients

Authors: Alexander D. Gilbert, Ivan G. Graham, Frances Y. Kuo, Robert Scheichl, Ian H. Sloan

Abstract: 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 b… ▽ More 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. △ Less

Submitted 17 May, 2019; v1 submitted 8 August, 2018; originally announced August 2018.

A Stein variational Newton method

Authors: Gianluca Detommaso, Tiangang Cui, Alessio Spantini, Youssef Marzouk, Robert Scheichl

Abstract: Stein variational gradient descent (SVGD) was recently proposed as a general purpose nonparametric variational inference algorithm [Liu & Wang, NIPS 2016]: it minimizes the Kullback-Leibler divergence between the target distribution and its approximation by implementing a form of functional gradient descent on a reproducing kernel Hilbert space. In this paper, we accelerate and generalize the SVGD… ▽ More Stein variational gradient descent (SVGD) was recently proposed as a general purpose nonparametric variational inference algorithm [Liu & Wang, NIPS 2016]: it minimizes the Kullback-Leibler divergence between the target distribution and its approximation by implementing a form of functional gradient descent on a reproducing kernel Hilbert space. In this paper, we accelerate and generalize the SVGD algorithm by including second-order information, thereby approximating a Newton-like iteration in function space. We also show how second-order information can lead to more effective choices of kernel. We observe significant computational gains over the original SVGD algorithm in multiple test cases. △ Less

Submitted 29 October, 2018; v1 submitted 8 June, 2018; originally announced June 2018.

Comments: 18 pages, 7 figures

Journal ref: NIPS 2018

Continuous Level Monte Carlo and Sample-Adaptive Model Hierarchies

Authors: Gianluca Detommaso, Tim Dodwell, Rob Scheichl

Abstract: In this paper, we present a generalisation of the Multilevel Monte Carlo (MLMC) method to a setting where the level parameter is a continuous variable. This Continuous Level Monte Carlo (CLMC) estimator provides a natural framework in PDE applications to adapt the model hierarchy to each sample. In addition, it can be made unbiased with respect to the expected value of the true quantity of interes… ▽ More In this paper, we present a generalisation of the Multilevel Monte Carlo (MLMC) method to a setting where the level parameter is a continuous variable. This Continuous Level Monte Carlo (CLMC) estimator provides a natural framework in PDE applications to adapt the model hierarchy to each sample. In addition, it can be made unbiased with respect to the expected value of the true quantity of interest provided the quantity of interest converges sufficiently fast. The practical implementation of the CLMC estimator is based on interpolating actual evaluations of the quantity of interest at a finite number of resolutions. As our new level parameter, we use the logarithm of a goal-oriented finite element error estimator for the accuracy of the quantity of interest. We prove the unbiasedness, as well as a complexity theorem that shows the same rate of complexity for CLMC as for MLMC. Finally, we provide some numerical evidence to support our theoretical results, by successfully testing CLMC on a standard PDE test problem. The numerical experiments demonstrate clear gains for sample-wise adaptive refinement strategies over uniform refinements. △ Less

Submitted 21 February, 2018; originally announced February 2018.

Comments: 22 pages, 4 figures

Circulant embedding with QMC -- analysis for elliptic PDE with lognormal coefficients

Authors: Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens, Rob Scheichl, Ian H. Sloan

Abstract: 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 t… ▽ More 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. △ Less

Submitted 2 April, 2018; v1 submitted 25 October, 2017; originally announced October 2017.

MSC Class: 60G10; 60G60; 65C05; 65C60; 35Q86; 65D32

Analysis of circulant embedding methods for sampling stationary random fields

Authors: Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens, Rob Scheichl, Ian H. Sloan

Abstract: 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 fie… ▽ More 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. △ Less

Submitted 20 March, 2018; v1 submitted 2 October, 2017; originally announced October 2017.

MSC Class: 60G10; 60G60; 65C05; 65C60

A hybrid Alternating Least Squares -- TT Cross algorithm for parametric PDEs

Authors: Sergey Dolgov, Robert Scheichl

Abstract: We consider the approximate solution of parametric PDEs using the low-rank Tensor Train (TT) decomposition. Such parametric PDEs arise for example in uncertainty quantification problems in engineering applications. We propose an algorithm that is a hybrid of the alternating least squares and the TT cross methods. It computes a TT approximation of the whole solution, which is beneficial when multip… ▽ More We consider the approximate solution of parametric PDEs using the low-rank Tensor Train (TT) decomposition. Such parametric PDEs arise for example in uncertainty quantification problems in engineering applications. We propose an algorithm that is a hybrid of the alternating least squares and the TT cross methods. It computes a TT approximation of the whole solution, which is beneficial when multiple quantities of interest are sought. This might be needed, for example, for the computation of the probability density function (PDF) via the maximum entropy method [Kavehrad and Joseph, IEEE Trans. Comm., 1986]. The new algorithm exploits and preserves the block diagonal structure of the discretized operator in stochastic collocation schemes. This disentangles computations of the spatial and parametric degrees of freedom in the TT representation. In particular, it only requires solving independent PDEs at a few parameter values, thus allowing the use of existing high performance PDE solvers. In our numerical experiments, we apply the new algorithm to the stochastic diffusion equation and compare it with preconditioned steepest descent in the TT format, as well as with (multilevel) quasi-Monte Carlo and dimension-adaptive sparse grids methods. For sufficiently smooth random fields the new approach is orders of magnitude faster. △ Less

Submitted 5 July, 2018; v1 submitted 14 July, 2017; originally announced July 2017.

Comments: 29 pages

MSC Class: 65F10; 65F30; 65N22; 65N30; 65N35

dune-composites -- A New Framework for High-Performance Finite Element Modelling of Laminates

Authors: Anne Reinarz, Tim Dodwell, Tim Fletcher, Linus Seelinger, Richard Butler, Robert Scheichl

Abstract: Finite element (FE) analysis has the potential to offset much of the expensive experimental testing currently required to certify aerospace laminates. However, large numbers of degrees of freedom are necessary to model entire aircraft components whilst accurately resolving micro-scale defects. The new module dune-composites, implemented within DUNE by the authors, provides a tool to efficiently so… ▽ More Finite element (FE) analysis has the potential to offset much of the expensive experimental testing currently required to certify aerospace laminates. However, large numbers of degrees of freedom are necessary to model entire aircraft components whilst accurately resolving micro-scale defects. The new module dune-composites, implemented within DUNE by the authors, provides a tool to efficiently solve large-scale problems using novel iterative solvers. The key innovation is a preconditioner that guarantees a constant number of iterations regardless of the problem size. Its robustness has been shown rigorously in Spillane et al. (Numer. Math. 126, 2014) for isotropic problems. For anisotropic problems in composites it is verified numerically for the first time in this paper. The parallel implementation in DUNE scales almost optimally over thousands of cores. To demonstrate this, we present an original numerical study, varying the shape of a localised wrinkle and the effect this has on the strength of a curved laminate. This requires a high-fidelity mesh containing at least four layers of quadratic elements across each ply and interface layer, underlining the need for dune-composites, which can achieve run times of just over 2 minutes on 2048 cores for realistic composites problems with 173 million degrees of freedom. △ Less

Submitted 13 July, 2017; originally announced July 2017.

Journal ref: Composite Structures, Volume 184, 2018, Pages 269-278

Modern Monte Carlo Variants for Uncertainty Quantification in Neutron Transport

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

Abstract: 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 include… ▽ More 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. △ Less

Submitted 17 October, 2017; v1 submitted 12 February, 2017; originally announced February 2017.

MSC Class: 65N06; 65N75; 65C05; 65C30; 65C50

Multilevel Monte Carlo and Improved Timestepping Methods in Atmospheric Dispersion Modelling

Authors: Grigoris Katsiolides, Eike H. Müller, Robert Scheichl, Tony Shardlow, Michael B. Giles, David J. Thomson

Abstract: A common way to simulate the transport and spread of pollutants in the atmosphere is via stochastic Lagrangian dispersion models. Mathematically, these models describe turbulent transport processes with stochastic differential equations (SDEs). The computational bottleneck is the Monte Carlo algorithm, which simulates the motion of a large number of model particles in a turbulent velocity field; f… ▽ More A common way to simulate the transport and spread of pollutants in the atmosphere is via stochastic Lagrangian dispersion models. Mathematically, these models describe turbulent transport processes with stochastic differential equations (SDEs). The computational bottleneck is the Monte Carlo algorithm, which simulates the motion of a large number of model particles in a turbulent velocity field; for each particle, a trajectory is calculated with a numerical timestepping method. Choosing an efficient numerical method is particularly important in operational emergency-response applications, such as tracking radioactive clouds from nuclear accidents or predicting the impact of volcanic ash clouds on international aviation, where accurate and timely predictions are essential. In this paper, we investigate the application of the Multilevel Monte Carlo (MLMC) method to simulate the propagation of particles in a representative one-dimensional dispersion scenario in the atmospheric boundary layer. MLMC can be shown to result in asymptotically superior computational complexity and reduced computational cost when compared to the Standard Monte Carlo (StMC) method, which is currently used in atmospheric dispersion modelling. To reduce the absolute cost of the method also in the non-asymptotic regime, it is equally important to choose the best possible numerical timestepping method on each level. To investigate this, we also compare the standard symplectic Euler method, which is used in many operational models, with two improved timestepping algorithms based on SDE splitting methods. △ Less

Submitted 20 September, 2017; v1 submitted 22 December, 2016; originally announced December 2016.

Comments: 30 pages, 13 figures, 3 tables, revised version, submitted to Journal of Computational Physics

MSC Class: 60H35; 65C05; 65C30; 76F40; 86A10 ACM Class: G.1.7; G.3; G.4; J.2

Scheduling massively parallel multigrid for multilevel Monte Carlo methods

Authors: Björn Gmeiner, Daniel Drzisga, Ulrich Ruede, Robert Scheichl, Barbara Wohlmuth

Abstract: The computational complexity of naive, sampling-based uncertainty quantification for 3D partial differential equations is extremely high. Multilevel approaches, such as multilevel Monte Carlo (MLMC), can reduce the complexity significantly, but to exploit them fully in a parallel environment, sophisticated scheduling strategies are needed. Often fast algorithms that are executed in parallel are es… ▽ More The computational complexity of naive, sampling-based uncertainty quantification for 3D partial differential equations is extremely high. Multilevel approaches, such as multilevel Monte Carlo (MLMC), can reduce the complexity significantly, but to exploit them fully in a parallel environment, sophisticated scheduling strategies are needed. Often fast algorithms that are executed in parallel are essential to compute fine level samples in 3D, whereas to compute individual coarse level samples only moderate numbers of processors can be employed efficiently. We make use of multiple instances of a parallel multigrid solver combined with advanced load balancing techniques. In particular, we optimize the concurrent execution across the three layers of the MLMC method: parallelization across levels, across samples, and across the spatial grid. The overall efficiency and performance of these methods will be analyzed. Here the scalability window of the multigrid solver is revealed as being essential, i.e., the property that the solution can be computed with a range of process numbers while maintaining good parallel efficiency. We evaluate the new scheduling strategies in a series of numerical tests, and conclude the paper demonstrating large 3D scaling experiments. △ Less

Submitted 12 July, 2016; originally announced July 2016.

MSC Class: G.1.8

Quasi-Monte Carlo and Multilevel Monte Carlo Methods for Computing Posterior Expectations in Elliptic Inverse Problems

Authors: R. Scheichl, A. M. Stuart, A. L. Teckentrup

Abstract: We are interested in computing the expectation of a functional of a PDE solution under a Bayesian posterior distribution. Using Bayes' rule, we reduce the problem to estimating the ratio of two related prior expectations. For a model elliptic problem, we provide a full convergence and complexity analysis of the ratio estimator in the case where Monte Carlo, quasi-Monte Carlo or multilevel Monte Ca… ▽ More We are interested in computing the expectation of a functional of a PDE solution under a Bayesian posterior distribution. Using Bayes' rule, we reduce the problem to estimating the ratio of two related prior expectations. For a model elliptic problem, we provide a full convergence and complexity analysis of the ratio estimator in the case where Monte Carlo, quasi-Monte Carlo or multilevel Monte Carlo methods are used as estimators for the two prior expectations. We show that the computational complexity of the ratio estimator to achieve a given accuracy is the same as the corresponding complexity of the individual estimators for the numerator and the denominator. We {also include numerical simulations, in the context of the model elliptic problem, which demonstrate the effectiveness of the approach. △ Less

Submitted 2 March, 2017; v1 submitted 15 February, 2016; originally announced February 2016.

Robust Numerical Upscaling of Elliptic Multiscale Problems at High Contrast

Authors: Daniel Peterseim, Robert Scheichl

Abstract: We present a new approach to the numerical upscaling for elliptic problems with rough diffusion coefficient at high contrast. It is based on the localizable orthogonal decomposition of $H^1$ into the image and the kernel of some novel stable quasi-interpolation operators with local $L^2$-approximation properties, independent of the contrast. We identify a set of sufficient assumptions on these qua… ▽ More We present a new approach to the numerical upscaling for elliptic problems with rough diffusion coefficient at high contrast. It is based on the localizable orthogonal decomposition of $H^1$ into the image and the kernel of some novel stable quasi-interpolation operators with local $L^2$-approximation properties, independent of the contrast. We identify a set of sufficient assumptions on these quasi-interpolation operators that guarantee in principle optimal convergence without pre-asymptotic effects for high-contrast coefficients. We then give an example of a suitable operator and establish the assumptions for a particular class of high-contrast coefficients. So far this is not possible without any pre-asymptotic effects, but the optimal convergence is independent of the contrast and the asymptotic range is largely improved over other discretisation schemes. The new framework is sufficiently flexible to allow also for other choices of quasi-interpolation operators and the potential for fully robust numerical upscaling at high contrast. △ Less

Submitted 25 January, 2016; originally announced January 2016.