Mark Ainsworth

Mark Ainsworth, Guosheng Fu

We present a new set of basis functions for H(curl)-conforming, H(div)-conforming, and L2 -conforming finite elements of arbitrary order on a tetrahedron. The basis functions are expressed in terms of Bernstein polynomials and augment the natural H1 -conforming Bernstein basis. The basis functions respect the differential operators, namely, the gradients of the high-order H1 -conforming Bernstein-Bezier basis functions form part of the H(curl)-conforming basis, and the curl of the high-order, non-gradients H(curl)-conforming basis functions form part of the H(div)-conforming basis, and the divergence of the high-order, non-curl H(div)-conforming basis functions form part of the L2-conforming basis. Procedures are given for the efficient computation of the mass and stiffness matrices with these basis functions without using quadrature rules for (piece-wise) constant coefficients on affine tetrahedra. Numerical results are presented to illustrate the use of the basis to approximate representative problems.

Mark Ainsworth, Johnny Guzmán, Francisco-Javier Sayas

The existence of uniformly bounded discrete extension operators is established for conforming Raviart-Thomas and Nédélec discretisations of $H(div)$ and $H(curl)$ on locally refined partitions of a polyhedral domain into tetrahedra.

Mark Ainsworth, Tomas Vejchodsky

A simple flux reconstruction for finite element solutions of reaction-diffusion problems is shown to yield fully computable upper bounds on the energy norm of error in an approximation of singularly perturbed reaction-diffusion problem. The flux reconstruction is based on simple, independent post-processing operations over patches of elements in conjunction with standard Raviart--Thomas vector fields and gives upper bounds even in cases where Galerkin orthogonality might be violated. If Galerkin orthogonality holds, we prove that the corresponding local error indicators are locally efficient and robust with respect to any mesh size and any size of the reaction coefficient, including the singularly perturbed limit.

Mark Ainsworth, Guosheng Fu

Composite basis functions for pyramidal elements on the spaces $H^1(Ω)$, $H(\mathrm{curl},Ω)$, $H(\mathrm{div},Ω)$ and $L^2(Ω)$ are presented. In particular, we construct the lowest-order composite pyramidal elements and show that they respect the de Rham diagram, i.e. we have an exact sequence and satisfy the commuting property. Moreover, the finite elements are fully compatible with the standard finite elements for the lowest-order Raviart-Thomas-Nédélec sequence on tetrahedral and hexahedral elements. That is to say, the new elements have the same degrees of freedom on the shared interface with the neighbouring hexahedral or tetrahedra elements, and the basis functions are conforming in the sense that they maintain the required level of continuity (full, tangential component, normal component, ...) across the interface. Furthermore, we study the approximation properties of the spaces as an initial partition consisting of tetrahedra, hexahedra and pyramid elements is successively subdivided and show that the spaces result in the same (optimal) order of approximation in terms of the mesh size $h$ as one would obtain using purely hexahedral or purely tetrahedral partitions.

Mark Ainsworth, Christian Glusa

The predicted reduced resiliency of next-generation high performance computers means that it will become necessary to take into account the effects of randomly occurring faults on numerical methods. Further, in the event of a hard fault occurring, a decision has to be made as to what remedial action should be taken in order to resume the execution of the algorithm. The action that is chosen can have a dramatic effect on the performance and characteristics of the scheme. Ideally, the resulting algorithm should be subjected to the same kind of mathematical analysis that was applied to the original, deterministic variant. The purpose of this work is to provide an analysis of the behaviour of the multigrid algorithm in the presence of faults. Multigrid is arguably the method of choice for the solution of large-scale linear algebra problems arising from discretization of partial differential equations and it is of considerable importance to anticipate its behaviour on an exascale machine. The analysis of resilience of algorithms is in its infancy and the current work is perhaps the first to provide a mathematical model for faults and analyse the behaviour of a state-of-the-art algorithm under the model. It is shown that the Two Grid Method fails to be resilient to faults. Attention is then turned to identifying the minimal necessary remedial action required to restore the rate of convergence to that enjoyed by the ideal fault-free method.

Mark Ainsworth, Guosheng Fu

The dispersive behavior of the recently proposed energy-conserving discontinuous Galerkin (DG) method by Fu and Shu [10] is analyzed and compared with the classical centered and upwinding DG schemes. It is shown that the new scheme gives a significant improvement over the classical centered and upwinding DG schemes in terms of dispersion error. Numerical results are presented to support the theoretical findings.

Mark Ainsworth, Justin Dong

We present a new approach to using neural networks to approximate the solutions of variational equations, based on the adaptive construction of a sequence of finite-dimensional subspaces whose basis functions are realizations of a sequence of neural networks. The finite-dimensional subspaces are then used to define a standard Galerkin approximation of the variational equation. This approach enjoys a number of advantages, including: the sequential nature of the algorithm offers a systematic approach to enhancing the accuracy of a given approximation; the sequential enhancements provide a useful indicator for the error that can be used as a criterion for terminating the sequential updates; the basic approach is largely oblivious to the nature of the partial differential equation under consideration; and, some basic theoretical results are presented regarding the convergence (or otherwise) of the method which are used to formulate basic guidelines for applying the method.

Mark Ainsworth, Yeonjong Shin

The ability of neural networks to provide `best in class' approximation across a wide range of applications is well-documented. Nevertheless, the powerful expressivity of neural networks comes to naught if one is unable to effectively train (choose) the parameters defining the network. In general, neural networks are trained by gradient descent type optimization methods, or a stochastic variant thereof. In practice, such methods result in the loss function decreases rapidly at the beginning of training but then, after a relatively small number of steps, significantly slow down. The loss may even appear to stagnate over the period of a large number of epochs, only to then suddenly start to decrease fast again for no apparent reason. This so-called plateau phenomenon manifests itself in many learning tasks. The present work aims to identify and quantify the root causes of plateau phenomenon. No assumptions are made on the number of neurons relative to the number of training data, and our results hold for both the lazy and adaptive regimes. The main findings are: plateaux correspond to periods during which activation patterns remain constant, where activation pattern refers to the number of data points that activate a given neuron; quantification of convergence of the gradient flow dynamics; and, characterization of stationary points in terms solutions of local least squares regression lines over subsets of the training data. Based on these conclusions, we propose a new iterative training method, the Active Neuron Least Squares (ANLS), characterised by the explicit adjustment of the activation pattern at each step, which is designed to enable a quick exit from a plateau. Illustrative numerical examples are included throughout.

Mark Ainsworth, Christian Glusa

A numerical scheme is presented for approximating fractional order Poisson problems in two and three dimensions. The scheme is based on reformulating the original problem posed over $Ω$ on the extruded domain $\mathcal{C}=Ω\times[0,\infty)$ following Caffarelli and Silvestre (2007). The resulting degenerate elliptic integer order PDE is then approximated using a hybrid FEM-spectral scheme. Finite elements are used in the direction parallel to the problem domain $Ω$, and an appropriate spectral method is used in the extruded direction. The spectral part of the scheme requires that we approximate the true eigenvalues of the integer order Laplacian over $Ω$. We derive an a priori error estimate which takes account of the error arising from using an approximation in place of the true eigenvalues. We further present a strategy for choosing approximations of the eigenvalues based on Weyl's law and finite element discretizations of the eigenvalue problem. The system of linear algebraic equations arising from the hybrid FEM-spectral scheme is decomposed into blocks which can be solved effectively using standard iterative solvers such as multigrid and conjugate gradient. Numerical examples in two and three dimensions show that the approach is quasi-optimal in terms of complexity.

Mark Ainsworth, Christian Glusa

We develop all of the components needed to construct an adaptive finite element code that can be used to approximate fractional partial differential equations, on non-trivial domains in $d\geq 1$ dimensions. Our main approach consists of taking tools that have been shown to be effective for adaptive boundary element methods and, where necessary, modifying them so that they can be applied to the fractional PDE case. Improved a priori error estimates are derived for the case of quasi-uniform meshes which are seen to deliver sub-optimal rates of convergence owing to the presence of singularities. Attention is then turned to the development of an a posteriori error estimate and error indicators which are suitable for driving an adaptive refinement procedure. We assume that the resulting refined meshes are locally quasi-uniform and develop efficient methods for the assembly of the resulting linear algebraic systems and their solution using iterative methods, including the multigrid method. The storage of the dense matrices along with efficient techniques for computing the dense matrix vector products needed for the iterative solution is also considered. The performance and efficiency of the resulting algorithm is illustrated for a variety of examples.

Mark Ainsworth, Christian Glusa

We explore the connection between fractional order partial differential equations in two or more spatial dimensions with boundary integral operators to develop techniques that enable one to efficiently tackle the integral fractional Laplacian. In particular, we develop techniques for the treatment of the dense stiffness matrix including the computation of the entries, the efficient assembly and storage of a sparse approximation and the efficient solution of the resulting equations. The main idea consists of generalising proven techniques for the treatment of boundary integral equations to general fractional orders. Importantly, the approximation does not make any strong assumptions on the shape of the underlying domain and does not rely on any special structure of the matrix that could be exploited by fast transforms. We demonstrate the flexibility and performance of this approach in a couple of two-dimensional numerical examples.

Mark Ainsworth, Guosheng Fu

We derive a posteriori error estimates for the hybridizable discontinuous Galerkin (HDG) methods, including both the primal and mixed formulations, for the approximation of a linear second-order elliptic problem on conforming simplicial meshes in two and three dimensions. We obtain fully computable, constant free, a posteriori error bounds on the broken energy seminorm and the HDG energy (semi)norm of the error. The estimators are also shown to provide local lower bounds for the HDG energy (semi)norm of the error up to a constant and a higher-order data oscillation term. For the primal HDG methods and mixed HDG methods with an appropriate choice of stabilization parameter, the estimators are also shown to provide a lower bound for the broken energy seminorm of the error up to a constant and a higher-order data oscillation term. Numerical examples are given illustrating the theoretical results.

Mark Ainsworth, Christian Glusa

Computing at the exascale level is expected to be affected by a significantly higher rate of faults, due to increased component counts as well as power considerations. Therefore, current day numerical algorithms need to be reexamined as to determine if they are fault resilient, and which critical operations need to be safeguarded in order to obtain performance that is close to the ideal fault-free method. In a previous paper, a framework for the analysis of random stationary linear iterations was presented and applied to the two grid method. The present work is concerned with the multigrid algorithm for the solution of linear systems of equations, which is widely used on high performance computing systems. It is shown that the Fault-Prone Multigrid Method is not resilient, unless the prolongation operation is protected. Strategies for fault detection and mitigation as well as protection of the prolongation operation are presented and tested, and a guideline for an optimal choice of parameters is devised.

Mark Ainsworth, Manuel A. Sánchez

The Bernstein-Bézier form of a polynomial is widely used in the fields of computer aided geometric design, spline approximation theory and, more recently, for high order finite element methods for the solution of partial differential equations. However, if one wishes to compute the classical Lagrange interpolant relative to the Bernstein basis, then the resulting Bernstein-Vandermonde matrix is found to be highly ill-conditioned. In the univariate case of degree $n$, Marco and Martinez showed that using Neville elimination to solve the system exploits the total positivity of the Bernstein basis and results in an $\mathcal{O}(n^2)$ complexity algorithm. Remarkable as it may be, the Marco-Martinez algorithm has some drawbacks: The derivation of the algorithm is quite technical; the interplay between the ideas of total positivity and Neville elimination are not part of the standard armoury of many non-specialists; and, the algorithm is strongly associated to the univariate case. The present work addresses these issues. An alternative algorithm for the inversion of the univariate Bernstein-Vandermonde system is presented that has: The same complexity as the Marco-Martinez algorithm and whose stability does not seem to be in any way inferior; a simple derivation using only the basic theory of Lagrange interpolation (at least in the univariate case); and, a natural generalisation to the multivariate case.

Mark Ainsworth, Tomáš Vejchodský

We present a fully computable a posteriori error estimator for piecewise linear finite element approximations of reaction-diffusion problems with mixed boundary conditions and piecewise constant reaction coefficient formulated in arbitrary dimension. The estimator provides a guaranteed upper bound on the energy norm of the error and it is robust for all values of the reaction coefficient, including the singularly perturbed case. The approach is based on robustly equilibrated boundary flux functions of Ainsworth and Oden (Wiley 2000) and on subsequent robust and explicit flux reconstruction. This paper simplifies and extends the applicability of the previous result of Ainsworth and Vejchodský (Numer. Math. 119 (2011) 219-243) in three aspects: (i) arbitrary dimension, (ii) mixed boundary conditions, and (iii) non-constant reaction coefficient. It is the first robust upper bound on the error with these properties. An auxiliary result that is of independent interest is the derivation of new explicit constants for two types of trace inequalities on simplices.