A Dynamically Adaptive Sparse Grid Method for Quasi-Optimal Interpolation of Multidimensional Analytic Functions

In this work we develop a dynamically adaptive sparse grids (SG) method for quasi-optimal interpolation of multidimensional analytic functions defined over a product of one dimensional bounded domains. The goal of such approach is to construct an interpolant in space that corresponds to the "best $M$-terms" based on sharp a priori estimate of polynomial coefficients. In the past, SG methods have been successful in achieving this, with a traditional construction that relies on the solution to a Knapsack problem: only the most profitable hierarchical surpluses are added to the SG. However, this approach requires additional sharp estimates related to the size of the analytic region and the norm of the interpolation operator, i.e., the Lebesgue constant. Instead, we present an iterative SG procedure that adaptively refines an estimate of the region and accounts for the effects of the Lebesgue constant. Our approach does not require any a priori knowledge of the analyticity or operator norm, is easily generalized to both affine and non-affine analytic functions, and can be applied to sparse grids build from one dimensional rules with arbitrary growth of the number of nodes. In several numerical examples, we utilize our dynamically adaptive SG to interpolate quantities of interest related to the solutions of parametrized elliptic and hyperbolic PDEs, and compare the performance of our quasi-optimal interpolant to several alternative SG schemes

An implementation of radiative transfer in the cosmological simulation code GADGET

We present a novel numerical implementation of radiative transfer in the cosmological smoothed particle hydrodynamics (SPH) simulation code {\small GADGET}. It is based on a fast, robust and photon-conserving integration scheme where the radiation transport problem is approximated in terms of moments of the transfer equation and by using a variable Eddington tensor as a closure relation, following the `OTVET'-suggestion of Gnedin & Abel. We derive a suitable anisotropic diffusion operator for use in the SPH discretization of the local photon transport, and we combine this with an implicit solver that guarantees robustness and photon conservation. This entails a matrix inversion problem of a huge, sparsely populated matrix that is distributed in memory in our parallel code. We solve this task iteratively with a conjugate gradient scheme. Finally, to model photon sink processes we consider ionisation and recombination processes of hydrogen, which is represented with a chemical network that is evolved with an implicit time integration scheme. We present several tests of our implementation, including single and multiple sources in static uniform density fields with and without temperature evolution, shadowing by a dense clump, and multiple sources in a static cosmological density field. All tests agree quite well with analytical computations or with predictions from other radiative transfer codes, except for shadowing. However, unlike most other radiative transfer codes presently in use for studying reionisation, our new method can be used on-the-fly during dynamical cosmological simulation, allowing simultaneous treatments of galaxy formation and the reionisation process of the Universe.Comment: 21 pages, 17 figures, published in MNRA

A Discrete Adapted Hierarchical Basis Solver For Radial Basis Function Interpolation

In this paper we develop a discrete Hierarchical Basis (HB) to efficiently solve the Radial Basis Function (RBF) interpolation problem with variable polynomial order. The HB forms an orthogonal set and is adapted to the kernel seed function and the placement of the interpolation nodes. Moreover, this basis is orthogonal to a set of polynomials up to a given order defined on the interpolating nodes. We are thus able to decouple the RBF interpolation problem for any order of the polynomial interpolation and solve it in two steps: (1) The polynomial orthogonal RBF interpolation problem is efficiently solved in the transformed HB basis with a GMRES iteration and a diagonal, or block SSOR preconditioner. (2) The residual is then projected onto an orthonormal polynomial basis. We apply our approach on several test cases to study its effectiveness, including an application to the Best Linear Unbiased Estimator regression problem

Adaptive Uncertainty Quantification for Stochastic Hyperbolic Conservation Laws

We propose a predictor-corrector adaptive method for the study of hyperbolic partial differential equations (PDEs) under uncertainty. Constructed around the framework of stochastic finite volume (SFV) methods, our approach circumvents sampling schemes or simulation ensembles while also preserving fundamental properties, in particular hyperbolicity of the resulting systems and conservation of the discrete solutions. Furthermore, we augment the existing SFV theory with a priori convergence results for statistical quantities, in particular push-forward densities, which we demonstrate through numerical experiments. By linking refinement indicators to regions of the physical and stochastic spaces, we drive anisotropic refinements of the discretizations, introducing new degrees of freedom (DoFs) where deemed profitable. To illustrate our proposed method, we consider a series of numerical examples for non-linear hyperbolic PDEs based on Burgers' and Euler's equations

Lowest order Virtual Element approximation of magnetostatic problems

We give here a simplified presentation of the lowest order Serendipity Virtual Element method, and show its use for the numerical solution of linear magneto-static problems in three dimensions. The method can be applied to very general decompositions of the computational domain (as is natural for Virtual Element Methods) and uses as unknowns the (constant) tangential component of the magnetic field $\mathbf{H}$ on each edge, and the vertex values of the Lagrange multiplier $p$ (used to enforce the solenoidality of the magnetic induction $\mathbf{B}=\mu\mathbf{H}$). In this respect the method can be seen as the natural generalization of the lowest order Edge Finite Element Method (the so-called "first kind N\'ed\'elec" elements) to polyhedra of almost arbitrary shape, and as we show on some numerical examples it exhibits very good accuracy (for being a lowest order element) and excellent robustness with respect to distortions

Comparison of some Reduced Representation Approximations

In the field of numerical approximation, specialists considering highly complex problems have recently proposed various ways to simplify their underlying problems. In this field, depending on the problem they were tackling and the community that are at work, different approaches have been developed with some success and have even gained some maturity, the applications can now be applied to information analysis or for numerical simulation of PDE's. At this point, a crossed analysis and effort for understanding the similarities and the differences between these approaches that found their starting points in different backgrounds is of interest. It is the purpose of this paper to contribute to this effort by comparing some constructive reduced representations of complex functions. We present here in full details the Adaptive Cross Approximation (ACA) and the Empirical Interpolation Method (EIM) together with other approaches that enter in the same category

Efficient Methods for Multidimensional Global Polynomial Approximation with Applications to Random PDEs

In this work, we consider several ways to overcome the challenges associated with polynomial approximation and integration of smooth functions depending on a large number of inputs. We are motivated by the problem of forward uncertainty quantification (UQ), whereby inputs to mathematical models are considered as random variables. With limited resources, finding more efficient and accurate ways to approximate the multidimensional solution to the UQ problem is of crucial importance, due to the “curse of dimensionality” and the cost of solving the underlying deterministic problem. The first way we overcome the complexity issue is by exploiting the structure of the approximation schemes used to solve the random partial differential equations (PDE), thereby significantly reducing the overall cost of the approximation. We do this first using multilevel approximations in the physical variables, and second by exploiting the hierarchy of nested sparse grids in the random parameter space. With these algorithmic advances, we provably decrease the complexity of collocation methods for solving random PDE problems. The second major theme in this work is the choice of efficient points for multidimensional interpolation and interpolatory quadrature. A major consideration in interpolation in multiple dimensions is the balance between stability, i.e., the Lebesgue constant of the interpolant, and the granularity of the approximation, e.g., the ability to choose an arbitrary number of interpolation points or to adaptively refine the grid. For these reasons, the Leja points are a popular choice for approximation on both bounded and unbounded domains. Mirroring the best-known results for interpolation on compact domains, we show that Leja points, defined for weighted interpolation on R, have a Lebesgue constant which grows subexponentially in the number of interpolation nodes. Regarding multidimensional quadratures, we show how certain new rules, generated from conformal mappings of classical interpolatory rules, can be used to increase the efficiency in approximating multidimensional integrals. Specifically, we show that the convergence rate for the novel mapped sparse grid interpolatory quadratures is improved by a factor that is exponential in the dimension of the underlying integral