Superconvergent interpolatory HDG methods for reaction diffusion equations I: An HDGk_{k} method

In our earlier work [8], we approximated solutions of a general class of scalar parabolic semilinear PDEs by an interpolatory hybridizable discontinuous Galerkin (Interpolatory HDG) method. This method reduces the computational cost compared to standard HDG since the HDG matrices are assembled once before the time integration. Interpolatory HDG also achieves optimal convergence rates; however, we did not observe superconvergence after an element-by-element postprocessing. In this work, we revisit the Interpolatory HDG method for reaction diffusion problems, and use the postprocessed approximate solution to evaluate the nonlinear term. We prove this simple change restores the superconvergence and keeps the computational advantages of the Interpolatory HDG method. We present numerical results to illustrate the convergence theory and the performance of the method

Interpolatory HDG Method for Parabolic Semilinear PDEs

We propose the interpolatory hybridizable discontinuous Galerkin (Interpolatory HDG) method for a class of scalar parabolic semilinear PDEs. The Interpolatory HDG method uses an interpolation procedure to efficiently and accurately approximate the nonlinear term. This procedure avoids the numerical quadrature typically required for the assembly of the global matrix at each iteration in each time step, which is a computationally costly component of the standard HDG method for nonlinear PDEs. Furthermore, the Interpolatory HDG interpolation procedure yields simple explicit expressions for the nonlinear term and Jacobian matrix, which leads to a simple unified implementation for a variety of nonlinear PDEs. For a globally-Lipschitz nonlinearity, we prove that the Interpolatory HDG method does not result in a reduction of the order of convergence. We display 2D and 3D numerical experiments to demonstrate the performance of the method

Analysis of \u3ci\u3ea posteriori\u3c/i\u3e error estimates of the discontinuous Galerkin method for nonlinear ordinary differential equations

We develop and analyze a new residual-based a posteriori error estimator for the discontinuous Galerkin (DG) method for nonlinear ordinary differential equations (ODEs). The a posteriori DG error estimator under investigation is computationally simple, efficient, and asymptotically exact. It is obtained by solving a local residual problem with no boundary condition on each element. We first prove that the DG solution exhibits an optimal O(hp+1) convergence rate in the L2-norm when p-degree piece-wise polynomials with p ≥1 are used. We further prove that the DG solution is O(h2p+1) superconvergent at the downwind points. We use these results to prove that the p-degree DG solution is O(hp+2) super close to a particular projection of the exact solution. This superconvergence result allows us to show that the true error can be divided into a significant part and a less significant part. The significant part of the discretization error for the DG solution is proportional to the (p +1)-degree right Radau polynomial and the less significant part converges at O(hp+2) rate in the L2-norm. Numerical experiments demonstrate that the theoretical rates are optimal. Based on the global superconvergent approximations, we construct asymptotically exact a posteriori error estimates and prove that they converge to the true errors in the L2-norm under mesh refinement. The order of convergence is proved to be p +2. Finally, we prove that the global effectivity index in the L2-norm converges to unity at O(h)rate. Several numerical examples are provided to illustrate the global superconvergence results and the convergence of the proposed estimator under mesh refinement. A local adaptive procedure that makes use of our local a posteriori error estimate is also presented

A uniformly accurate finite elements method for singular perturbation problems

We consider piecewise polynomial finite elements method for a singular perturbation problem. The finite elements method of Griffiths for a problem with non-constant coefficients was adapted by introducing piecewise polynomial approximation. We generate the tridiagonal difference schemes which are second order accurate in uniform norm

HDGlab: An Open-Source Implementation of the Hybridisable Discontinuous Galerkin Method in MATLAB

This paper presents HDGlab, an open source MATLAB implementation of the hybridisable discontinuous Galerkin (HDG) method. The main goal is to provide a detailed description of both the HDG method for elliptic problems and its implementation available in HDGlab. Ultimately, this is expected to make this relatively new advanced discretisation method more accessible to the computational engineering community. HDGlab presents some features not available in other implementations of the HDG method that can be found in the free domain. First, it implements high-order polynomial shape functions up to degree nine, with both equally-spaced and Fekete nodal distributions. Second, it supports curved isoparametric simplicial elements in two and three dimensions. Third, it supports non-uniform degree polynomial approximations and it provides a flexible structure to devise degree adaptivity strategies. Finally, an interface with the open-source high-order mesh generator Gmsh is provided to facilitate its application to practical engineering problems

Non-axisymmetric instabilities in self-gravitating tori around black holes, and solving Einstein constraints with superconvergent finite element methods

This thesis contains results on two related projects. In the first project, we explore non-axisymmetric instabilities in general relativistic accretion disks around black holes. Such disks are created as transient structures in several astrophysical scenarios, including mergers of compact objects and core collapse of massive stars. These disks are suggested for the role of cenral engines of gamma-ray bursts. We address the stability of these objects against the runaway and non-axisymmetric instabilities in the three-dimensional hydrodynamical fully general relativistic treatment. We explore three slender and moderately slender disk models with varying disk-to-black hole mass ratio. None of the models that we consider develop the runaway instability during the time span of the simulations, despite large radial axisymmetric oscillations, induced in the disks by the initial data construction procedure. All models develop unstable non-axisymmetric modes on a dynamical timescale. In simulations with dynamical general relativistic treatment, we observe two distinct types of instabilities: the Papaloizou-Pringle instability and the so-called Intermediate instability. The development of the nonaxisymmetric mode with azimuthal number m=1 is enhanced by the outspiraling motion of the black hole. The overall picture of the unstable modes in our disk models is similar to the Newtonian case. In the second project, we experiment with solving the Einstein constraint equations using finite elements on semistructured triangulations of multiblock grids. We illustrate our approach with a simple example of Brill wave initial data, with the constraints reducing to a single linear elliptic equation for the conformal factor $\psi$. We use quadratic Lagrange elements on semi-structured simplicial meshes, obtained by triangulation of multi-block grids. In the case of uniform refinement the scheme is superconvergent at most mesh vertices, due to local symmetry of the finite element basis with respect to local spatial inversions. We show that in the superconvergent case subsequent unstructured mesh refinements do not improve the quality of our initial data. As proof of concept that this approach is feasible for generating multi-block initial data in three dimensions, we evolve the constructed initial data using a high order finite-differencing multi-block approach and extract gravitational waves from the numerical solution

Layer-adapted meshes for convection-diffusion problems

This is a book on numerical methods for singular perturbation problems - in particular stationary convection-dominated convection-diffusion problems. More precisely it is devoted to the construction and analysis of layer-adapted meshes underlying these numerical methods. An early important contribution towards the optimization of numerical methods by means of special meshes was made by N.S. Bakhvalov in 1969. His paper spawned a lively discussion in the literature with a number of further meshes being proposed and applied to various singular perturbation problems. However, in the mid 1980s this development stalled, but was enlivend again by G.I. Shishkin's proposal of piecewise- equidistant meshes in the early 1990s. Because of their very simple structure they are often much easier to analyse than other meshes, although they give numerical approximations that are inferior to solutions on competing meshes. Shishkin meshes for numerous problems and numerical methods have been studied since and they are still very much in vogue. With this contribution we try to counter this development and lay the emphasis on more general meshes that - apart from performing better than piecewise-uniform meshes - provide a much deeper insight in the course of their analysis. In this monograph a classification and a survey are given of layer-adapted meshes for convection-diffusion problems. It tries to give a comprehensive review of state-of-the art techniques used in the convergence analysis for various numerical methods: finite differences, finite elements and finite volumes. While for finite difference schemes applied to one-dimensional problems a rather complete convergence theory for arbitrary meshes is developed, the theory is more fragmentary for other methods and problems and still requires the restriction to certain classes of meshes

A Multigrid Method for the Efficient Numerical Solution of Optimization Problems Constrained by Partial Differential Equations

We study the minimization of a quadratic functional subject to constraints given by a linear or semilinear elliptic partial differential equation with distributed control. Further, pointwise inequality constraints on the control are accounted for. In the linear-quadratic case, the discretized optimality conditions yield a large, sparse, and indefinite system with saddle point structure. One main contribution of this thesis consists in devising a coupled multigrid solver which avoids full constraint elimination. To this end, we define a smoothing iteration incorporating elements from constraint preconditioning. A local mode analysis shows that for discrete optimality systems, we can expect smoothing rates close to those obtained with respect to the underlying constraint PDE. Our numerical experiments include problems with constraints where standard pointwise smoothing is known to fail for the underlying PDE. In particular, we consider anisotropic diffusion and convection-diffusion problems. The framework of our method allows to include line smoothers or ILU-factorizations, which are suitable for such problems. In all cases, numerical experiments show that convergence rates do not depend on the mesh size of the finest level and discrete optimality systems can be solved with a small multiple of the computational cost which is required to solve the underlying constraint PDE. Employing the full multigrid approach, the computational cost is proportional to the number of unknowns on the finest grid level. We discuss the role of the regularization parameter in the cost functional and show that the convergence rates are robust with respect to both the fine grid mesh size and the regularization parameter under a mild restriction on the next to coarsest mesh size. Incorporating spectral filtering for the reduced Hessian in the control smoothing step allows us to weaken the mesh size restriction. As a result, problems with near-vanishing regularization parameter can be treated efficiently with a negligible amount of additional computational work. For fine discretizations, robust convergence is obtained with rates which are independent of the regularization parameter, the coarsest mesh size, and the number of levels. In order to treat linear-quadratic problems with pointwise inequality constraints on the control, the multigrid approach is modified to solve subproblems generated by a primal-dual active set strategy (PDAS). Numerical experiments demonstrate the high efficiency of this approach due to mesh-independent convergence of both the outer PDAS method and the inner multigrid solver. The PDAS-multigrid method is incorporated in the sequential quadratic programming (SQP) framework. Inexact Newton techniques further enhance the computational efficiency. Globalization is implemented with a line search based on the augmented Lagrangian merit function. Numerical experiments highlight the efficiency of the resulting SQP-multigrid approach. In all cases, locally superlinear convergence of the SQP method is observed. In combination with the mesh-independent convergence rate of the inner solver, a solution method with optimal efficiency is obtained

Error analysis of the Galerkin FEM in L 2 -based norms for problems with layers: On the importance, conception and realization of balancing

In the present thesis it is shown that the most natural choice for a norm for the analysis of the Galerkin FEM, namely the energy norm, fails to capture the boundary layer functions arising in certain reaction-diffusion problems. In view of a formal Definition such reaction-diffusion problems are not singularly perturbed with respect to the energy norm. This observation raises two questions: 1. Does the Galerkin finite element method on standard meshes yield satisfactory approximations for the reaction-diffusion problem with respect to the energy norm? 2. Is it possible to strengthen the energy norm in such a way that the boundary layers are captured and that it can be reconciled with a robust finite element method, i.e.~robust with respect to this strong norm? In Chapter 2 we answer the first question. We show that the Galerkin finite element approximation converges uniformly in the energy norm to the solution of the reaction-diffusion problem on standard shape regular meshes. These results are completely new in two dimensions and are confirmed by numerical experiments. We also study certain convection-diffusion problems with characterisitc layers in which some layers are not well represented in the energy norm. These theoretical findings, validated by numerical experiments, have interesting implications for adaptive methods. Moreover, they lead to a re-evaluation of other results and methods in the literature. In 2011 Lin and Stynes were the first to devise a method for a reaction-diffusion problem posed in the unit square allowing for uniform a priori error estimates in an adequate so-called balanced norm. Thus, the aforementioned second question is answered in the affirmative. Obtaining a non-standard weak formulation by testing also with derivatives of the test function is the key idea which is related to the H^1-Galerkin methods developed in the early 70s. Unfortunately, this direct approach requires excessive smoothness of the finite element space considered. Lin and Stynes circumvent this problem by rewriting their problem into a first order system and applying a mixed method. Now the norm captures the layers. Therefore, they need to be resolved by some layer-adapted mesh. Lin and Stynes obtain optimal error estimates with respect to the balanced norm on Shishkin meshes. However, their method is unable to preserve the symmetry of the problem and they rely on the Raviart-Thomas element for H^div-conformity. In Chapter 4 of the thesis a new continuous interior penalty (CIP) method is present, embracing the approach of Lin and Stynes in the context of a broken Sobolev space. The resulting method induces a balanced norm in which uniform error estimates are proven. In contrast to the mixed method the CIP method uses standard Q_2-elements on the Shishkin meshes. Both methods feature improved stability properties in comparison with the Galerkin FEM. Nevertheless, the latter also yields approximations which can be shown to converge to the true solution in a balanced norm uniformly with respect to diffusion parameter. Again, numerical experiments are conducted that agree with the theoretical findings. In every finite element analysis the approximation error comes into play, eventually. If one seeks to prove any of the results mentioned on an anisotropic family of Shishkin meshes, one will need to take advantage of the different element sizes close to the boundary. While these are ideally suited to reflect the solution behavior, the error analysis is more involved and depends on anisotropic interpolation error estimates. In Chapter 3 the beautiful theory of Apel and Dobrowolski is extended in order to obtain anisotropic interpolation error estimates for macro-element interpolation. This also sheds light on fundamental construction principles for such operators. The thesis introduces a non-standard finite element space that consists of biquadratic C^1-finite elements on macro-elements over tensor product grids, which can be viewed as a rectangular version of the C^1-Powell-Sabin element. As an application of the general theory developed, several interpolation operators mapping into this FE space are analyzed. The insight gained can also be used to prove anisotropic error estimates for the interpolation operator induced by the well-known C^1-Bogner-Fox-Schmidt element. A special modification of Scott-Zhang type and a certain anisotropic interpolation operator are also discussed in detail. The results of this chapter are used to approximate the solution to a recation-diffusion-problem on a Shishkin mesh that features highly anisotropic elements. The obtained approximation features continuous normal derivatives across certain edges of the mesh, enabling the analysis of the aforementioned CIP method.:Notation 1 Introduction 2 Galerkin FEM error estimation in weak norms 2.1 Reaction-diffusion problems 2.2 A convection-diffusion problem with weak characteristic layers and a Neumann outflow condition 2.3 A mesh that resolves only part of the exponential layer and neglects the weaker characteristic layers 2.3.1 Weakly imposed characteristic boundary conditions 2.4 Numerical experiments 2.4.1 A reaction-diffusion problem with boundary layers 2.4.2 A reaction-diffusion problem with an interior layer 2.4.3 A convection-diffusion problem with characteristic layers and a Neumann outflow condition 2.4.4 A mesh that resolves only part of the exponential layer and neglects the weaker characteristic layers 3 Macro-interpolation on tensor product meshes 3.1 Introduction 3.2 Univariate C1-P2 macro-element interpolation 3.3 C1-Q2 macro-element interpolation on tensor product meshes 3.4 A theory on anisotropic macro-element interpolation 3.5 C1 macro-interpolation on anisotropic tensor product meshes 3.5.1 A reduced macro-element interpolation operator 3.5.2 The full C1-Q2 interpolation operator 3.5.3 A C1-Q2 macro-element quasi-interpolation operator of Scott-Zhang type on tensor product meshes 3.5.4 Summary: anisotropic C1 (quasi-)interpolation error estimates 3.6 An anisotropic macro-element of tensor product type 3.7 Application of macro-element interpolation on a tensor product Shishkin mesh 4 Balanced norm results for reaction-diffusion 4.1 The balanced finite element method of Lin and Stynes 4.2 A C0 interior penalty method 4.3 Galerkin finite element method 4.3.1 L2-norm error bounds and supercloseness 4.3.2 Maximum-norm error bounds 4.4 Numerical verification 4.5 Further developments and summary Reference