A piecewise linear finite element discretization of the diffusion equation for arbitrary polyhedral grids

We develop a piecewise linear (PWL) Galerkin finite element spatial discretization for the multi-dimensional radiation diffusion equation. It uses recently introduced piecewise linear weight and basis functions in the finite element approximation and it can be applied on arbitrary polygonal (2D) or polyhedral (3D) grids. We first demonstrate some analytical properties of the PWL method and perform a simple mode analysis to compare the PWL method with Palmer’s vertex-centered finite-volume method and with a bilinear continuous finite element method. We then show that this new PWL method gives solutions comparable to those from Palmer’s. However, since the PWL method produces a symmetric positive-definite coefficient matrix, it should be substantially more computationally efficient than Palmer’s method, which produces an asymmetric matrix. We conclude that the Galerkin PWL method is an attractive option for solving diffusion equations on unstructured grids.     

EXPLICIT PREDICTOR-MULTICORRECTOR TIME DISCONTINUOUS GALERKIN METHODS FOR NON-LINEAR DYNAMICS

Explicit predictor-multicorrector time discontinuous Galerkin (TDG) methods developed for linear structural dynamics are formulated and implemented in a form suitable for arbitrary non-linear analysis of structural dynamics problems. The formulation is intended to inherit the accuracy properties of the exact parent implicit TDG methods. To this end, suitable predictors and correctors are designed to achieve third order accuracy, large stability limits and controllable numerical dissipation by means of an algorithmic parameter. As the study of a general non-linear case is rather complex, the analysis of the convergence properties of the resulting algorithms are restricted to conservative Duffing oscillators, for which closed-form solutions are available. It is shown that the main properties of the underlying parent scheme can be retained. Finally, results of representative numerical simulations relevant to Duffing oscillators and to a stiff spring pendulum discretized with finite elements illustrate the performance of the numerical schemes and confirm the analytical estimates.     

Explicit one-step time discretizations for discontinuous Galerkin and finite volume schemes based on local predictors

We consider a family of explicit one-step time discretizations for finite volume and discontinuous Galerkin schemes, which is based on a predictor-corrector formulation. The predictor remains local taking into account the time evolution of the data only within the grid cell. Based on a space–time Taylor expansion, this idea is already inherent in the MUSCL finite volume scheme to get second order accuracy in time and was generalized in the context of higher order ENO finite volume schemes. We interpret the space–time Taylor expansion used in this approach as a local predictor and conclude that other space–time approximate solutions of the local Cauchy problem in the grid cell may be applied. Three possibilities are considered in this paper: (1) the classical space–time Taylor expansion, in which time derivatives are obtained from known space-derivatives by the Cauchy–Kovalewsky procedure; (2) a local continuous extension Runge–Kutta scheme and (3) a local space–time Galerkin predictor with a version suitable for stiff source terms. The advantage of the predictor–corrector formulation is that the time evolution is done in one step which establishes optimal locality during the whole time step. This time discretization scheme can be used within all schemes which are based on a piecewise continuous approximation as finite volume schemes, discontinuous Galerkin schemes or the recently proposed reconstructed discontinuous Galerkin or P_NP_M schemes. The implementation of these approaches is described, advantages and disadvantages of different predictors are discussed and numerical results are shown.     

Central discontinuous Galerkin methods for ideal MHD equations with the exactly divergence-free magnetic field

In this paper, central discontinuous Galerkin methods are developed for solving ideal magnetohydrodynamic (MHD) equations. The methods are based on the original central discontinuous Galerkin methods designed for hyperbolic conservation laws on overlapping meshes, and use different discretization for magnetic induction equations. The resulting schemes carry many features of standard central discontinuous Galerkin methods such as high order accuracy and being free of exact or approximate Riemann solvers. And more importantly, the numerical magnetic field is exactly divergence-free. Such property, desired in reliable simulations of MHD equations, is achieved by first approximating the normal component of the magnetic field through discretizing induction equations on the mesh skeleton, namely, the element interfaces. And then it is followed by an element-by-element divergence-free reconstruction with the matching accuracy. Numerical examples are presented to demonstrate the high order accuracy and the robustness of the schemes.     

Finite Difference/Element Method for a Two-Dimensional Modified Fractional Diffusion Equation

We present the finite difference/element method for a two-dimensional modified fractional diffusion equation. The analysis is carried out first for the time semi-discrete scheme, and then for the full discrete scheme. The time discretization is based on the $L1$-approximation for the fractional derivative terms and the second-order backward differentiation formula for the classical first order derivative term. We use finite element method for the spatial approximation in full discrete scheme. We show that both the semi-discrete and full discrete schemes are unconditionally stable and convergent. Moreover, the optimal convergence rate is obtained. Finally, some numerical examples are tested in the case of one and two space dimensions and the numerical results confirm our theoretical analysis.     

Krylov implicit integration factor methods for spatial discretization on high dimensional unstructured meshes: Application to discontinuous Galerkin methods

Integration factor methods are a class of “exactly linear part” time discretization methods. In [Q. Nie, Y.-T. Zhang, R. Zhao, Efficient semi-implicit schemes for stiff systems, Journal of Computational Physics, 214 (2006) 521–537], a class of efficient implicit integration factor (IIF) methods were developed for solving systems with both stiff linear and nonlinear terms, arising from spatial discretization of time-dependent partial differential equations (PDEs) with linear high order terms and stiff lower order nonlinear terms. The tremendous challenge in applying IIF temporal discretization for PDEs on high spatial dimensions is how to evaluate the matrix exponential operator efficiently. For spatial discretization on unstructured meshes to solve PDEs on complex geometrical domains, how to efficiently apply the IIF temporal discretization was open. In this paper, we solve this problem by applying the Krylov subspace approximations to the matrix exponential operator. Then we apply this novel time discretization technique to discontinuous Galerkin (DG) methods on unstructured meshes for solving reaction–diffusion equations. Numerical examples are shown to demonstrate the accuracy, efficiency and robustness of the method in resolving the stiffness of the DG spatial operator for reaction–diffusion PDEs. Application of the method to a mathematical model in pattern formation during zebrafish embryo development shall be shown.     

Local DG method using WENO type limiters for convection–diffusion problems

The local discontinuous Galerkin (LDG) method is a spatial discretization procedure for convection–diffusion equations, which employs useful features from high resolution finite volume schemes, such as the exact or approximate Riemann solvers serving as numerical fluxes and limiters, which is termed as Runge–Kutta LDG (RKLDG) when TVD Runge–Kutta method is applied for time discretization. It has the advantage of flexibility in handling complicated geometry, h-p adaptivity, and efficiency of parallel implementation and has been used successfully in many applications. However, the limiters used to control spurious oscillations in the presence of strong shocks are less robust than the strategies of essentially non-oscillatory (ENO) and weighted ENO (WENO) finite volume and finite difference methods. In this paper, we investigated RKLDG methods with WENO and Hermite WENO (HWENO) limiters for solving convection–diffusion equations on unstructured meshes, with the goal of obtaining a robust and high order limiting procedure to simultaneously obtain uniform high order accuracy and sharp, non-oscillatory shock transition. Numerical results are provided to illustrate the behavior of these procedures.     

Stability and accuracy of spatial approximations for wave equation tidal models

Amplitude and phase characteristics for numerical approximations to the shallow water wave equation are obtained for linear and quadratic finite elements, for finite difference approximations, for non-constant bathemetry, and for uneven node spacing. Stability is shown to require non-zero friction as well as satisfaction of a Courant constraint. Lumping is shown to reduce the Courant constraint for stability while higher order and quadratic finite element approximations require a more restrictive constraint than their second order and linear finite element counterparts. The amplitude of the propagation factor for stable schemes and propagating waves is seen to be independent of the Courant number and type of numerical approximation. Although the higher order and quadratic schemes provide better propagation of the low and moderate frequency waves, the highest frequency waves (2Δx) are better propagated by low order numerical methods.     

A class of hybrid DG/FV methods for conservation laws I: Basic formulation and one-dimensional systems

By comparing the discontinuous Galerkin (DG) and the finite volume (FV) methods, a concept of ‘static reconstruction’ and ‘dynamic reconstruction’ is introduced for high-order numerical methods. Based on the new concept, a class of hybrid DG/FV schemes is presented for one-dimensional conservation law using a ‘hybrid reconstruction’ approach. In the hybrid DG/FV schemes, the lower-order derivatives of a piecewise polynomial solution are computed locally in a cell by the DG method based on Taylor basis functions (called as ‘dynamic reconstruction’), while the higher-order derivatives are re-constructed by the ‘static reconstruction’ of the FV method, using the known lower-order derivatives in the cell itself and its adjacent neighboring cells. The hybrid DG/FV methods can greatly reduce CPU time and memory required by the traditional DG methods with the same order of accuracy on the same mesh, and they can be extended directly to unstructured and hybrid grids in two and three dimensions similar to the DG and/or FV methods. The hybrid DG/FV methods are applied to one-dimensional conservation law, including linear and non-linear scalar equation and Euler equations. In order to capture the strong shock waves without spurious oscillations, a simple shock detection approach is developed to mark ‘trouble cells’, and a moment limiter is adopted for higher-order schemes. The numerical results demonstrate the accuracy, and the super-convergence property is shown for the third-order hybrid DG/FV schemes. In addition, by analyzing the eigenvalues of the semi-discretized system in one dimension, we discuss the spectral properties of the hybrid DG/FV schemes to explain the super-convergence phenomenon.     

Comparison of solvers for the generalized Riemann problem for hyperbolic systems with source terms

We compare four different approximate solvers for the generalized Riemann problem (GRP) for non-linear systems of hyperbolic equations with source terms. The GRP is a special Cauchy problem for a hyperbolic system with source terms whose initial condition is piecewise smooth. We briefly review the four solvers currently available and carry out a systematic assessment of these in terms of accuracy and computational efficiency. These solvers are the building block for constructing high-order numerical schemes of the ADER type for solving the general initial-boundary value problem for inhomogeneous systems in multiple space dimensions, in the frameworks of finite volume and discontinuous Galerkin finite element methods.     

High-order upwind residual distribution schemes on isoparametric curved elements

Residual distribution schemes on curved geometries are discussed in the context of higher order spatial discretization for hyperbolic conservation laws. The discrete solution, defined by a Finite Element space based on triangular Lagrangian P_k elements, is globally continuous. A natural sub-triangulation of these elements allows to reuse the simple distribution schemes previously developed for linear P₁ triangles. The paper introduces curved elements with piecewise quadratic and cubic approximation of the boundaries of the domain, using standard sub- or isoparametric transformation. Numerical results for the Euler equations confirm the predicted order of accuracy, showing the importance of a higher order approximation of the geometry.     

A discontinuous Galerkin method for inviscid low Mach number flows

In this work we extend the high-order discontinuous Galerkin (DG) finite element method to inviscid low Mach number flows. The method here presented is designed to improve the accuracy and efficiency of the solution at low Mach numbers using both explicit and implicit schemes for the temporal discretization of the compressible Euler equations. The algorithm is based on a classical preconditioning technique that in general entails modifying both the instationary term of the governing equations and the dissipative term of the numerical flux function (full preconditioning approach). In the paper we show that full preconditioning is beneficial for explicit time integration while the implicit scheme turns out to be efficient and accurate using just the modified numerical flux function. Thus the implicit scheme could also be used for time accurate computations. The performance of the method is demonstrated by solving an inviscid flow past a NACA0012 airfoil at different low Mach numbers using various degrees of polynomial approximations. Computations with and without preconditioning are performed on different grid topologies to analyze the influence of the spatial discretization on the accuracy of the DG solutions at low Mach numbers.     

Bad behavior of Godunov mixed methods for strongly anisotropic advection–dispersion equations

We study the performance of Godunov mixed methods, which combine a mixed-hybrid finite element solver and a Godunov-like shock-capturing solver, for the numerical treatment of the advection–dispersion equation with strong anisotropic tensor coefficients. It turns out that a mesh locking phenomenon may cause ill-conditioning and reduce the accuracy of the numerical approximation especially on coarse meshes. This problem may be partially alleviated by substituting the mixed-hybrid finite element solver used in the discretization of the dispersive (diffusive) term with a linear Galerkin finite element solver, which does not display such a strong ill conditioning. To illustrate the different mechanisms that come into play, we investigate the spectral properties of such numerical discretizations when applied to a strongly anisotropic diffusive term on a small regular mesh. A thorough comparison of the stiffness matrix eigenvalues reveals that the accuracy loss of the Godunov mixed method is a structural feature of the mixed-hybrid method. In fact, the varied response of the two methods is due to the different way the smallest and largest eigenvalues of the dispersion (diffusion) tensor influence the diagonal and off-diagonal terms of the final stiffness matrix. One and two dimensional test cases support our findings.     

A class of hybrid DG/FV methods for conservation laws II: Two-dimensional cases

By comparing the discontinuous Galerkin (DG) methods, the k-exact finite volume (FV) methods and the lift collocation penalty (LCP) methods, a concept of ‘static reconstruction’ and ‘dynamic reconstruction’ was introduced for higher-order numerical methods in our previous work. Based on this concept, a class of hybrid DG/FV methods was presented for one-dimensional conservation law using a ‘hybrid reconstruction’ approach. In the hybrid DG/FV schemes, the lower-order derivatives of the piecewise polynomial are computed locally in a cell by the traditional DG method (called as ‘dynamic reconstruction’), while the higher-order derivatives are re-constructed by the ‘static reconstruction’ of the FV method, using the known lower-order derivatives in the cell itself and in its adjacent face neighboring cells. In this follow-up paper, the hybrid DG/FV schemes are extended onto two-dimensional unstructured and hybrid grids. The two-dimensional linear and non-linear scalar conservation law and Euler equations are considered. Some typical cases are tested to demonstrate the performance of the hybrid DG/FV method, and the numerical results show that they can reduce the CPU time and memory requirement greatly than the traditional DG method with the same order of accuracy in the same mesh.     

High order conservative Lagrangian schemes with Lax–Wendroff type time discretization for the compressible Euler equations

In this paper, we explore the Lax–Wendroff (LW) type time discretization as an alternative procedure to the high order Runge–Kutta time discretization adopted for the high order essentially non-oscillatory (ENO) Lagrangian schemes developed in 3 and 5. The LW time discretization is based on a Taylor expansion in time, coupled with a local Cauchy–Kowalewski procedure to utilize the partial differential equation (PDE) repeatedly to convert all time derivatives to spatial derivatives, and then to discretize these spatial derivatives based on high order ENO reconstruction. Extensive numerical examples are presented, for both the second-order spatial discretization using quadrilateral meshes [3] and third-order spatial discretization using curvilinear meshes [5]. Comparing with the Runge–Kutta time discretization procedure, an advantage of the LW time discretization is the apparent saving in computational cost and memory requirement, at least for the two-dimensional Euler equations that we have used in the numerical tests.     

Finite Element $θ$-Schemes for the Acoustic Wave Equation

In this paper, we investigate the stability and convergence of a family of implicit finite difference schemes in time and Galerkin finite element methods in space for the numerical solution of the acoustic wave equation. The schemes cover the classical explicit second-order leapfrog scheme and the fourth-order accurate scheme in time obtained by the modified equation method. We derive general stability conditions for the family of implicit schemes covering some well-known CFL conditions. Optimal error estimates are obtained. For sufficiently smooth solutions, we demonstrate that the maximal error in the $L^2$-norm error over a finite time interval converges optimally as $\mathcal{O}(h^{p+1}+∆t^s)$, where $p$ denotes the polynomial degree, $s$=2 or 4, $h$ the mesh size, and $∆t$ the time step.     

On maximum-principle-satisfying high order schemes for scalar conservation laws

We construct uniformly high order accurate schemes satisfying a strict maximum principle for scalar conservation laws. A general framework (for arbitrary order of accuracy) is established to construct a limiter for finite volume schemes (e.g. essentially non-oscillatory (ENO) or weighted ENO (WENO) schemes) or discontinuous Galerkin (DG) method with first order Euler forward time discretization solving one-dimensional scalar conservation laws. Strong stability preserving (SSP) high order time discretizations will keep the maximum principle. It is straightforward to extend the method to two and higher dimensions on rectangular meshes. We also show that the same limiter can preserve the maximum principle for DG or finite volume schemes solving two-dimensional incompressible Euler equations in the vorticity stream-function formulation, or any passive convection equation with an incompressible velocity field. Numerical tests for both the WENO finite volume scheme and the DG method are reported.     

A stable high-order Spectral Difference method for hyperbolic conservation laws on triangular elements

Numerical schemes using piecewise polynomial approximation are very popular for high order discretization of conservation laws. While the most widely used numerical scheme under this paradigm appears to be the Discontinuous Galerkin method, the Spectral Difference scheme has often been found attractive as well, because of its simplicity of formulation and implementation. However, recently it has been shown that the scheme is not linearly stable on triangles. In this paper we present an alternate formulation of the scheme, featuring a new flux interpolation technique using Raviart–Thomas spaces, which proves stable under a similar linear analysis in which the standard scheme failed. We demonstrate viability of the concept by showing linear stability both in the semi-discrete sense and for time stepping schemes of the SSP Runge–Kutta type. Furthermore, we present convergence studies, as well as case studies in compressible flow simulation using the Euler equations.     

Explicit Runge–Kutta residual distribution schemes for time dependent problems: Second order case

In this paper, we construct spatially consistent explicit second order discretizations for time dependent hyperbolic problems, starting from a given residual distribution (RD) discrete approximation of the steady operator. We review the existing knowledge on consistent RD mass matrices and highlight the relations between different definitions. We then introduce our explicit approach which is based on three main ingredients: first recast the RD discretization as a stabilized Galerkin scheme, then use a shifted time discretization in the stabilization operator, and lastly apply high order mass lumping on the Galerkin component of the discretization. The discussion is particularly relevant for schemes of the residual distribution type 18 and 3 which we will use for all our numerical experiments. However, similar ideas can be used in the context of residual-based finite volume discretizations such as the ones proposed in 14 and 12. The schemes are tested on a wide variety of classical problems confirming the theoretical expectations.     

Finite element method for modeling radiative transfer in semitransparent graded index cylindrical medium

Both Galerkin finite element method (GFEM) and least squares finite element method (LSFEM) are developed and their performances are compared for solving the radiative transfer equation of graded index medium in cylindrical coordinate system (RTEGC). The angular redistribution term of the RTEGC is discretized by finite difference approach and after angular discretization the RTEGC is formulated into a discrete-ordinates form, which is then discretized based on Galerkin or least squares finite element approach. To overcome the RTEGC-led numerical singularity at the origin of cylindrical coordinate system, a pole condition is proposed as a special mathematical boundary condition. Compared with the GFEM, the LSFEM has very good numerical properties and can effectively mitigate the nonphysical oscillation appeared in the GFEM solutions. Various problems of both axisymmetry and nonaxisymmetry, and with medium of uniform refractive index distribution or graded refractive index distribution are tested. The results show that both the finite element approaches have good accuracy to predict the radiative heat transfer in semitransparent graded index cylindrical medium, while the LSFEM has better numerical stability.