A survey of Trefftz methods for the Helmholtz equation

Trefftz methods are finite element-type schemes whose test and trial functions are (locally) solutions of the targeted differential equation. They are particularly popular for time-harmonic wave problems, as their trial spaces contain oscillating basis functions and may achieve better approximation properties than classical piecewise-polynomial spaces. We review the construction and properties of several Trefftz variational formulations developed for the Helmholtz equation, including least squares, discontinuous Galerkin, ultra weak variational formulation, variational theory of complex rays and wave based methods. The most common discrete Trefftz spaces used for this equation employ generalised harmonic polynomials (circular and spherical waves), plane and evanescent waves, fundamental solutions and multipoles as basis functions; we describe theoretical and computational aspects of these spaces, focusing in particular on their approximation properties. One of the most promising, but not yet well developed, features of Trefftz methods is the use of adaptivity in the choice of the propagation directions for the basis functions. The main difficulties encountered in the implementation are the assembly and the ill-conditioning of linear systems, we briefly survey some strategies that have been proposed to cope with these problems.Comment: 41 pages, 2 figures, to appear as a chapter in Springer Lecture Notes in Computational Science and Engineering. Differences from v1: added a few sentences in Sections 2.1, 2.2.2 and 2.3.1; inserted small correction

Development and Optimization of Non-Hydrostatic Models for Water Waves and Fluid-Vegetation Interaction

The primary objective of this study is twofold: 1) to develop an efficient and accurate non-hydrostatic wave model for fully dispersive highly nonlinear waves, and 2) to investigate the interaction between waves and submerged flexible vegetation using a fully coupled wave-vegetation model. This research consists of three parts. Firstly, an analytical dispersion relationship is derived for waves simulated by models utilizing Keller-box scheme and central differencing for vertical discretization. The phase speed can be expressed as a rational polynomial function of the dimensionless water depth, $kh$, and the layer distribution in water column becomes an optimizable parameter in this function. For a given tolerance dispersion error, the range of $kh$ is extended and the layer thicknesses are optimally selected. The derived theoretical dispersion relationship is tested with linear and nonlinear standing waves generated by an Euler model. The optimization method is applicable to other non-hydrostatic models for water waves. Secondly, an efficient and accurate approach is developed to solve Euler equations for fully dispersive and highly nonlinear water waves. Discontinuous Galerkin, finite difference, and spectral element formulations are used for horizontal discretization, vertical discretization, and the Poisson equation, respectively. The Keller-box scheme is adopted for its capability of resolving frequency dispersion accurately with low vertical resolution (two or three layers). A three-stage optimal Strong Stability-Preserving Runge-Kutta (SSP-RK) scheme is employed for time integration. Thirdly, a fully coupled wave-vegetation model for simulating the interaction between water waves and submerged flexible plants is presented. The complete governing equation for vegetation motion is solved with a high-order finite element method and an implicit time differencing scheme. The vegetation model is fully coupled with a wave model to explore the relationship between displacement of water particle and plant stem, as well as the effect of vegetation flexibility on wave attenuation. This vegetation deformation model can be coupled with other wave models to simulate wave-vegetation interactions

Bilinear Immersed Finite Elements for Interface Problems

In this dissertation we discuss bilinear immersed finite elements (IFE) for solving interface problems. The related research works can be categorized into three aspects: (1) the construction of the bilinear immersed finite element spaces; (2) numerical methods based on these IFE spaces for solving interface problems; and (3) the corresponding error analysis. All of these together form a solid foundation for the bilinear IFEs. The research on immersed finite elements is motivated by many real world applications, in which a simulation domain is often formed by several materials separated from each other by curves or surfaces while a mesh independent of interface instead of a body-fitting mesh is preferred. The bilinear IFE spaces are nonconforming finite element spaces and the mesh can be independent of interface. The error estimates for the interpolation of a Sobolev function in a bilinear IFE space indicate that this space has the usual approximation capability expected from bilinear polynomials, which is O(h2) in L2 norm and O(h) in H1 norm. Then the immersed spaces are applied in Galerkin, finite volume element (FVE) and discontinuous Galerkin (DG) methods for solving interface problems. Numerical examples show that these methods based on the bilinear IFE spaces have the same optimal convergence rates as those based on the standard bilinear finite element for solutions with certain smoothness. For the symmetric selective immersed discontinuous Galerkin method based on bilinear IFE, we have established its optimal convergence rate. For the Galerkin method based on bilinear IFE, we have also established its convergence. One of the important advantages of the discontinuous Galerkin method is its flexibility for both p and h mesh refinement. Because IFEs can use a mesh independent of interface, such as a structured mesh, the combination of a DG method and IFEs allows a flexible adaptive mesh independent of interface to be used for solving interface problems. That is, a mesh independent of interface can be refined wherever needed, such as around the interface and the singular source. We also develop an efficient selective immersed discontinuous Galerkin method. It uses the sophisticated discontinuous Galerkin formulation only around the locations needed, but uses the simpler Galerkin formulation everywhere else. This selective formulation leads to an algebraic system with far less unknowns than the immersed DG method without scarifying the accuracy; hence it is far more efficient than the conventional discontinuous Galerkin formulations

Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018

This open access book features a selection of high-quality papers from the presentations at the International Conference on Spectral and High-Order Methods 2018, offering an overview of the depth and breadth of the activities within this important research area. The carefully reviewed papers provide a snapshot of the state of the art, while the extensive bibliography helps initiate new research directions

Reconstruction based error detection for robust approximation of partial differential equations

In this work we present a framework for the construction of robust a posteriori estimates for classes of finite difference schemes. We are motivated by the relative lack of such frameworks compared to existing ones for other numerical discretisation methods, such as finite elements and finite volumes. The framework we propose is based on the use of reconstructions, which are obtained by post-processing the finite difference solution. The post-processed object is a key ingredient in obtaining a posteriori bounds using the relevant stability framework of the problem. The resulting bounds are fully computable and allow us to establish a posteriori error control over the problem at hand. In the first part of the thesis we motivate and investigate the behaviour of our framework using model ODE, elliptic and hyperbolic problems. We use our framework to obtain reconstructions which are used to compute a posteriori error estimates. We validate the numerical behaviour of these estimates using solutions of varying regularity. In the second part of the thesis we focus on hyperbolic conservation laws in one spatial dimension and we deal with scalar problems as well as systems. Hyperbolic conservation laws are widely used in the modelling of physical phenomena. The numerical modelling of conservation laws, which arises due to the frequent lack of explicit solutions, is challenging, largely due to the complex behaviour these problems exhibit, such as shock formation even with smooth initial conditions. In this setting, we present a framework which is applicable to general non-linear conservation laws. We investigate its numerical behaviour and showcase our results by using popular finite difference discretisations for a range of problems. We demonstrate that the the framework can produce optimal estimates, capable of tracking features of interest and act as refinement/coarsening indicators

On entropy satisfying and maximum-principle-satisfying high-order methods for Fokker-Planck equations

Computation of Fokker-Planck equations with satisfying long time behavior is important in many applications. In this thesis, we design, analyze and implement entropy satisfying and maximum-principle-satisfying high-order numerical methods to solve the Fokker-Planck equation of the finitely extensible nonlinear elastic (FENE) dumbbell model for polymers, subject to homogeneous fluids, and the reaction-diffusion-advection equation arising in the evolution of biased dispersal of population dynamics. The design of each method is guided to satisfy three main properties, consisting of the nonnegativity principle, the mass conservation and the preservation of nonzero steady states. The relative entropy and the maximum principle are two powerful tools used to evaluate our methods, for instance, the steady state preservation can be ensured if the method is either entropy satisfying or maximum principle satisfying in the sense that the ratio of the solution to the equilibrium will stay in the same range as indicated by the initial data. These schemes are constructed in several steps, including reformulation of the Fokker-Planck equation into its nonlogarithmic Landau form, spacial discretization by discontinuous Galerken (DG) methods and some Runge-Kutta time discretization. The special form of numerical fluxes motivated by those introduced in [H. Liu and J. Yan, Commun. Comput. Phys. 8(3), 2010, 541-564] is essential to incorporate desired properties into each scheme through choices of flux parameters. In this thesis, we have obtained the following results. 1. For the Fokker-Planck equation of the FENE model, we propose an entropy satisfying conservative method which preserves all the three desired properties at both semidiscrete and discrete levels. This method is shown to be entropy satisfying in the sense that these schemes satisfy discrete entropy inequalities for both the physical entropy and the quadratic entropy. These ensure that the computed solution is a probability density, and the schemes are entropy stable and preserve the equilibrium solutions. 2. We further develop an entropy satisfying DG method of arbitrary high order. Both semidiscrete and fully discrete methods are shown to satisfy two desired properties: mass conservation and entropy satisfying for the quadratic entropy, therefore preserving the equilibrium solutions. A positive numerical approximation is obtained with the same accuracy as the numerical solution through a reconstruction at the final time. For both the finite volume scheme and the DG scheme we also prove the convergence of numerical solutions to the equilibrium solution as time tends to infinity. One- and two-dimensional numerical results are provided to demonstrate the good qualities of these schemes and effects of some canonical homogeneous flows. 3. We develop up to third-order accurate DG methods satisfying a strict maximum principle for a class of linear Fokker-Planck equations. A procedure is established to identify an effective test set in each computational cell to ensure the desired bounds of numerical averages during time evolution. This is achievable by properly choosing flux parameters and a positive decomposition of weighted cell averages. Based on this result, a scaling limiter for the DG method with Euler forward time discretization is proposed to solve both one- and multidimensional Fokker-Planck equations. As a consequence, the present scheme preserves steady states and provides a satisfying long time behavior. Numerical tests for the DG method are reported, with applications to polymer models with both Hookean and FENE potentials. 4. For Fokker-Planck equations with reaction such as the biased dispersal model in population dynamics, we develop entropy/energy stable finite difference schemes. For the numerical method to capture the long-time pattern of persistence or extinction, we use the relative entropy when the resource potential is logarithmic and explore the usual energy for other resource potentials. The present schemes are shown to satisfy three important properties of the continuous model for the population density: (i) positivity preserving, (ii) equilibrium preserving and (iii) entropy or energy satisfying. These ensure that our schemes provide a satisfying long-time behavior, thus revealing the desired dispersal pattern. Moreover, we present several numerical results which confirm the second-order accuracy for various resource potentials and underline the efficiency to preserve the large time asymptotic

Numerical Analysis of Flux Reconstruction

High-order methods have become of increasing interest in recent years in computational physics. This is in part due to their perceived ability to, in some cases, reduce the computational overhead of complex problems through both an efficient use of computational resources and a reduction in the required degrees of freedom. One such high-order method in particular – Flux Reconstruction – is the focus of this thesis. This body of work relies and expands on the theoretical methods that are used to understand the behaviour of numerical methods – particularly related to their real-world application to industrial problems. The thesis begins by challenging some of the existing dogma surrounding computational fluid dynamics by evaluating the performance of high-order flux reconstruction. First, the use of the primitive variables as an intermediary step in the construction of flux terms is investigated. It is found that reducing the order of the flux function by using the conserved rather than primitive variables has a substantial impact on the resolution of the method. Critically, this is supported by a theoretical analysis, which shows that this mechanism of error generation becomes increasing important to consider as the order of accuracy increases. Next, the analysis of Flux Reconstruction was extended by analytically and numerically exploring the impact of higher dimensionality and grid deformation. It is found that both expanding and contracting grids – essential components of real-world domain decomposition – can cause dispersion overshoot in two dimensions, but that FR appears to suffer less that comparable Finite Difference approaches. Fully discrete analysis is then used to show that, depending on the correction function, small perturbations in incidence angle can cause large changes in group velocity. The same analysis is also used to theoretically demonstrate that Discontinuous Galerkin suffers less from dispersion error than Huynh’s FR scheme – a phenomenon that has previously been observed experimentally, but not explained theoretically. This thesis concludes with the presentation of a robust theoretical underpinning for determining stable correction functions for FR. Three new families of correction functions are presented, and their properties extensively explored. An important theoretical finding is introduced – that stable correction functions are not defined uniquely be a norm. As a result, a generalised approach is presented, which is able to recover all previously defined correction functions, but in some instances via a different norm to their original derivation. This new super-family of correction functions shows considerable promise in increasing temporal stability limits, reducing dispersion when fully discretised, and increasing the rate of convergence. Taken altogether, this thesis represents a considerable advance in the theoretical characterisation and understanding of a numerical method – one which, it has been shown, has enormous potential for forming the heart of future computational physics codes

Numerical analysis of a fluid droplet subject to acoustic waves

Efficient and rigorous acoustic solvers that enable high frequency sweep application over a wide range of frequencies are of great interest due to their practical importance in many engineering, physical problems or life science research that involve acoustic radiation, such as engine noise analysis, acoustic simulation in micro-fluidics and the design of lab device, etc. There is room for reduction of cost on experimental systems that can be investigated and optimised through numerical modelling of physical processes on the micro-scale level. The major difficulty that arises is the inconsistency of materials, time scales and fast oscillation nature of the solution that leads to unstable results for conventional numerical methods. However, analytical solutions are infeasible for large problems with complex geometries and sophisticated boundary conditions. Hence, the vital need for efficient solvers. In this research the development of computational methods for acoustic application is presented. The proposed method is applied to the study of propagating waves in particular to simulate acoustic phenomena in micro-droplet actuated by leaky Surface Acoustic Waves on a lithium niobate (LiNbO3) substrate. Explicitly, we introduce a new computational method for the analysis of fluids subjected to high frequency mechanical forcing. Here we solve the Helmholtz equation in the frequency domain, applying higher order Lobatto hierarchical finite element approximation in H1 space, where both pressure field and geometry are independently approximated with arbitrary and heterogeneous polynomial order. Meanwhile, a time dependent acoustic solver with arbitrary input signals is also proposed and implemented. The development of extended computational methods for the solution of the Helmholtz equation with polychromatic waves is presented, where Fourier transformation is applied to switch the incident wave and solution space from the frequency domain to the temporal domain. Consequently, the implementation and convergence rate of the numerical methods are demonstrated with benchmark problems. The numerical method is an extension of the conventional higher order finite element method and as such it relies on the definition of basis functions. In this work we implement a set of basis functions using integrated Legendre polynomials (Lobatto polynomial). Two type of basis functions are presented and compared. Therefore, the significant improvements in efficiency is demonstrated using a Lobatto hierarchical basis compared with a Legendre type basis. Moreover, a novel error estimation and automatic adaptivity scheme is outlined based on an existing a priori error estimator. The accuracy and efficiency of the proposed object oriented (predefined error level) a priori error estimator is validated through numerical assessments on a three-dimensional spherical problem and compared with uniformly h and p adaptivities. The simple and generic features of the proposed scheme allow fast frequency sweeps with low computational cost for multiple frequencies acoustic application. The current finite element approach is executed in parallel with pre-partitioned domain, which guarantees the optimal computational speed with minimal computational effort for large problems. Overall, the benefits of using the proposed acoustic solver is explained in detail. Finally, we illustrate the model's performance using an example of a micro-droplet actuated by a surface acoustic wave (SAW), which has vast applications in micro-fluidics and micro-rheology at high frequency. Conclusions are drawn, and future directions are pointed out. The proposed finite element technology is implemented in the University of Glasgow in-house open-source finite element parallel computational code, MoFEM (Mesh Oriented Finite Element Method). All algorithms and examples are publicly available for download and testing