Tensor Product Elements Research Articles

The Spectral Difference Raviart–Thomas Method for Two and Three-Dimensional Elements and Its Connection with the Flux Reconstruction Formulation

The purpose of this work is to describe in detail the development of the spectral difference Raviart–Thomas (SDRT) formulation for two and three-dimensional tensor-product elements and simplexes. Through the process, the authors establish the equivalence between the SDRT method and the flux reconstruction (FR) approach under the assumption of the linearity of the flux and the mesh uniformity. Such a connection allows building a new family of FR schemes for two and three-dimensional simplexes and also to recover the well-known FR-SD method with tensor-product elements. In addition, a thorough analysis of the numerical dissipation and dispersion of both aforementioned schemes and the nodal discontinuous Galerkin FR (FR-DG) method with two and three-dimensional elements is proposed through the use of the combined-mode Fourier approach. SDRT is shown to possess an enhanced temporal linear stability in comparison to FR-DG. On the contrary, SDRT displays larger dissipation and dispersion errors with respect to FR-DG. Finally, the study is concluded with a set of numerical experiments, the linear advection-diffusion problem, the Isentropic Euler Vortex, and the Taylor-Green Vortex (TGV). The latter test case shows that SDRT schemes present a non-linear unstable behavior with simplex elements and certain polynomial degrees. For the sake of completeness, the matrix form of the SDRT method is developed and the computational performance of SDRT with respect to FR schemes is evaluated using GPU architectures.

Anisotropic mesh generation and adaptation for quads using the Lp-CVT method

This paper presents a framework for anisotropic unstructured all-quad mesh adaptation. The technique is suitable for use with a diverse class of numerical methods, including the high-order Discontinuous Galerkin solver presented here. Benefiting from the flexibility and efficiency of these methods, improvements in the solution accuracy are targeted by modifying the discretization of the domain. Previous works have tackled this problem for the case of simplex elements using continuous mesh and error models. Our method makes use of novel extensions to allow similar techniques to be applied for the case of tensor product elements. Key aspects of the work will discuss the modified continuous mesh representation and error minimization procedure. Additionally, these ideas have been coupled with mesh generation based on both existing and experimental techniques. As a result, this work will also consider the implementation of an energy minimizing, node redistribution method, that takes place in an anisotropic spatially varying high-order p-norm. Finally, the technique will also be applied to goal-oriented adaption using the adjoint-based error estimation. Overall, the work outlines the adaptive procedure tailored specially to all-quad meshing.

Bringing Trimmed Serendipity Methods to Computational Practice in Firedrake

We present an implementation of the trimmed serendipity finite element family, using the open-source finite element package Firedrake. The new elements can be used seamlessly within the software suite for problems requiring H 1 , H (curl), or H (div)-conforming elements on meshes of squares or cubes. To test how well trimmed serendipity elements perform in comparison to traditional tensor product elements, we perform a sequence of numerical experiments including the primal Poisson, mixed Poisson, and Maxwell cavity eigenvalue problems. Overall, we find that the trimmed serendipity elements converge, as expected, at the same rate as the respective tensor product elements, while being able to offer significant savings in the time or memory required to solve certain problems.

Hyperbolic diffusion in flux reconstruction: Optimisation through kernel fusion within tensor-product elements

Novel methods are presented in this initial study for the fusion of GPU kernels in the artificial compressibility method (ACM), using tensor product elements with constant Jacobians and flux reconstruction. This is made possible through the hyperbolisation of the diffusion terms, which eliminates the expensive algorithmic steps needed to form the viscous stresses. Two fusion approaches are presented, which offer differing levels of parallelism. This is found to be necessary for the change in workload as the order of accuracy of the elements is increased. Several further optimisations of these approaches are demonstrated, including a generation time memory manager which maximises resource usage. The fused kernels are able to achieve 3-4 times speedup, which compares favourably with a theoretical maximum speedup of 4. In three dimensional test cases, the generated fused kernels are found to reduce total runtime by ∼25%, and, when compared to the standard ACM formulation, simulations demonstrate that a speedup of 2.3 times can be achieved.

Mortar-based Entropy-Stable Discontinuous Galerkin Methods on Non-conforming Quadrilateral and Hexahedral Meshes

High-order entropy-stable discontinuous Galerkin (DG) methods for nonlinear conservation laws reproduce a discrete entropy inequality by combining entropy conservative finite volume fluxes with summation-by-parts (SBP) discretization matrices. In the DG context, on tensor product (quadrilateral and hexahedral) elements, SBP matrices are typically constructed by collocating at Lobatto quadrature points. Recent work has extended the construction of entropy-stable DG schemes to collocation at more accurate Gauss quadrature points (Chan et al. in SIAM J Sci Comput 41(5):A2938–A2966, 2019) . In this work, we extend entropy-stable Gauss collocation schemes to non-conforming meshes. Entropy-stable DG schemes require computing entropy conservative numerical fluxes between volume and surface quadrature nodes. On conforming tensor product meshes where volume and surface nodes are aligned, flux evaluations are required only between “lines” of nodes. However, on non-conforming meshes, volume and surface nodes are no longer aligned, resulting in a larger number of flux evaluations. We reduce this expense by introducing an entropy-stable mortar-based treatment of non-conforming interfaces via a face-local correction term, and provide necessary conditions for high-order accuracy. Numerical experiments for the compressible Euler equations in two and three dimensions confirm the stability and accuracy of this approach.

Dispersion analysis of displacement-based and TDNNS mixed finite elements for thin-walled elastodynamics

We compare several lowest-order finite element approximations to the problem of elastodynamics of thin-walled structures by means of dispersion analysis, which relates the parameter frequency-times-thickness (fd) and the wave speed. We restrict to analytical theory of harmonic front-crested waves that freely propagate in an infinite plate. Our study is formulated as a quasi-periodic eigenvalue problem on a single tensor–product element, which is eventually layered in the thickness direction. In the first part of the paper it is observed that the displacement-based finite elements align with the theory provided there are sufficiently many layers. In the second part we present novel anisotropic hexahedral tangential-displacement and normal–normal-stress continuous (TDNNS) mixed finite elements for Hellinger–Reissner formulation of elastodynamics. It turns out that one layer of such elements is sufficient for fd up to 2000[kHzmm]. Nevertheless, due to a large amount of TDNNS degrees of freedom the computational complexity is only comparable to the multi-layer displacement-based element. This is not the case at low frequencies, where TDNNS is by far more efficient since it allows for rough anisotropic discretizations, contrary to the displacement-based elements that suffer from the shear locking effect.

Optimized geometrical metrics satisfying free-stream preservation

Computational fluid dynamics and aerodynamics, which complement more expensive empirical approaches, are critical for developing aerospace vehicles. During the past three decades, computational aerodynamics capability has improved remarkably, following advances in computer hardware and algorithm development. However, for complex applications, the demands on computational fluid dynamics continue to increase in a quest to gain a few percent improvements in accuracy. Herein, we numerically demonstrate, in the context of tensor-product discretizations on hexahedral elements, that computing the metric terms with an optimization-based approach leads to a solution whose accuracy is overall on par and often better than the one obtained using the widely adopted Thomas and Lombard metric terms computation (Geometric conservation law and its application to flow computations on moving grids, AIAA Journal, 1979). We show the efficacy of the proposed technique in the context of low and high-order accurate nonlinearly stable (entropy stable) schemes on distorted, high-order tensor product elements, considering smooth three-dimensional inviscid and viscous compressible test cases for which an analytical solution is known. The methodology, originally developed by Crean et al. (2018) in the context of triangular/tetrahedral grids, is not limited to tensor-product cells and it can be applied to other cell-based diagonal-norm summation-by-parts discretizations, including spectral differences, discontinuous Galerkin finite elements, and flux reconstruction schemes.

Hybrid multigrid methods for high-order discontinuous Galerkin discretizations

The present work develops hybrid multigrid methods for high-order discontinuous Galerkin discretizations of elliptic problems, which are, for example, a key ingredient of incompressible flow solvers in the field of computational fluid dynamics. Fast matrix-free operator evaluation on tensor product elements is used to devise a computationally efficient PDE solver. The multigrid hierarchy exploits all possibilities of geometric, polynomial, and algebraic coarsening, targeting engineering applications on complex geometries. Additionally, a transfer from discontinuous to continuous function spaces is performed within the multigrid hierarchy. This does not only further reduce the problem size of the coarse-grid problem, but also leads to a discretization most suitable for state-of-the-art algebraic multigrid methods applied as coarse-grid solver. The relevant design choices regarding the selection of optimal multigrid coarsening strategies among the various possibilities are discussed with the metric of computational costs as the driving force for algorithmic selections. We find that a transfer to a continuous function space at highest polynomial degree (or on the finest mesh), followed by polynomial and geometric coarsening, shows the best overall performance. The success of this particular multigrid strategy is due to a significant reduction in iteration counts as compared to a transfer from discontinuous to continuous function spaces at lowest polynomial degree (or on the coarsest mesh). The coarsening strategy with transfer to a continuous function space on the finest level leads to a multigrid algorithm that is robust with respect to the penalty parameter of the symmetric interior penalty method. Detailed numerical investigations are conducted for a series of examples ranging from academic test cases to more complex, practically relevant geometries. Performance comparisons to state-of-the-art methods from the literature demonstrate the versatility and computational efficiency of the proposed multigrid algorithms.

Efficient Matrix-Free High-Order Finite Element Evaluation for Simplicial Elements

With the gap between processor clock speeds and memory bandwidth speeds continuing to increase, the use of arithmetically intense schemes, such as high-order finite element methods, continues to be of considerable interest. In particular, the use of matrix-free formulations of finite element operators for tensor-product elements of quadrilaterals in two dimensions and hexahedra in three dimensions, in combination with single-instruction multiple-data instruction sets, is a well-studied topic at present for the efficient implicit solution of elliptic equations. However, a considerable limiting factor for this approach is the use of meshes comprising of only quadrilaterals or hexahedra, the creation of which is still an open problem within the mesh generation community. In this article, we study the efficiency of high-order finite element operators for the Helmholtz equation with a focus on extending this approach to unstructured meshes of triangles, tetrahedra, and prismatic elements using the spectral/$hp$ element method and corresponding tensor-product bases for these element types. We show that although performance is naturally degraded when going from hexahedra to these simplicial elements, efficient implementations can still be obtained that are capable of attaining 50% through 70% floating point operations of the peak of processors with both AVX2 and AVX512 instruction sets.

Efficient Entropy Stable Gauss Collocation Methods

The construction of high order entropy stable collocation schemes on quadrilateral and hexahedral elements has relied on the use of Gauss-Legendre-Lobatto collocation points and their equivalence with summation-by-parts (SBP) finite difference operators. In this work, we show how to efficiently generalize the construction of semi-discretely entropy stable schemes on tensor product elements to Gauss points and generalized SBP operators. Numerical experiments suggest that the use of Gauss points significantly improves accuracy on curved meshes.

On the exponent of exponential convergence of p-version FEM spaces

We study the exponent of the exponential rate of convergence in terms of the number of degrees of freedom for various non-standard p-version finite element spaces employing reduced cardinality basis. More specifically, we show that serendipity finite element methods and discontinuous Galerkin finite element methods with total degree mathcal {P}_{p} basis have a faster exponential convergence with respect to the number of degrees of freedom than their counterparts employing the tensor product mathcal {Q}_{p} basis for quadrilateral/hexahedral elements, for piecewise analytic problems under p-refinement. The above results are proven by using a new p-optimal error bound for the L2-orthogonal projection onto the total degree mathcal {P}_{p} basis, and for the H1-projection onto the serendipity finite element space over tensor product elements with dimension d ≥ 2. These new p-optimal error bounds lead to a larger exponent of the exponential rate of convergence with respect to the number of degrees of freedom. Moreover, these results show that part of the basis functions in mathcal {Q}_{p} basis plays no roles in achieving the hp-optimal error bound in the Sobolev space. The sharpness of theoretical results is also verified by a series of numerical examples.

Numerical studies of serendipity and tensor product elements for eigenvalue problems

While the use of finite element methods for the numerical approximation of eigenvalues is a well-studied problem, the use of serendipity elements for this purpose has received little attention in the literature. We show by numerical experiments that serendipity elements, which are defined on a square reference geometry, can attain the same order of accuracy as their tensor product counterparts while using dramatically fewer degrees of freedom. In some cases, the serendipity method uses only half as many basis functions as the tensor product method while still producing the same numerical approximation of an eigenvalue. To encourage the further use and study of serendipity elements, we provide a table of serendipity basis functions for low order cases and a Mathematica file that can be used to generate the basis functions for higher order cases.

An Analysis of Stability of the Flux Reconstruction Formulation on Quadrilateral Elements for the Linear Advection–Diffusion Equation

The Flux Reconstruction (FR) approach to high-order methods is a flexible and robust framework that has proven to be a promising alternative to the traditional Discontinuous Galerkin (DG) schemes on parallel architectures like Graphical Processing Units (GPUs) since it pairs exceptionally well with explicit time-stepping methods. The FR formulation was originally proposed by Huynh (AIAA Pap 2007-4079:1---42, 2007). Vincent et al. (J Sci Comput 47(1):50---72, 2011) later developed a single parameter family of correction functions which provide energy stable schemes under this formulation in 1D. These schemes, known as Vincent---Castonguay---Jameson---Huynh (VCJH) schemes, offer control over properties like stability, dispersion and dissipation through the variation of the VCJH parameter. Classical schemes like nodal-DG and Spectral Difference (SD) can also be recovered under this formulation. Following the development of the FR approach in 1D, Castonguay et al. (J Sci Comput 51(1):224---256, 2012) and Williams et al. (J Comput Phys 250:53---76, 2013) and Williams and Jameson (J Sci Comput 59(3):721---759, 2014) developed correction functions that give rise to energy stable FR formulations for triangles and tetrahedra. For the case of tensor product elements like quadrilaterals and hexahedra however, a simple extension of the 1D FR approach utilizing the 1D VCJH correction functions was possible and has been adopted by several authors (Castonguay in High-order energy stable flux reconstruction schemes for fluid flow simulations on unstructured grids, 2012; Witherden et al. in Comput Fluids 120:173---186, 2015; Comput Phys Commun 185(11):3028---3040, 2014). But whether such an extension of the 1D approach to tensor product elements is stable remained an open question. A direct extension of the 1D stability analysis fails due to certain key difficulties which necessitate the formulation of a norm different from the one utilized for stability analysis in 1D and on simplex elements. We have recently overcome these issues and shown that the VCJH schemes are stable for the linear advection equation on Cartesian meshes for any non-negative value of the VCJH parameter (Sheshadri and Jameson in J Sci Comput 67(2):769---790, 2016; J Sci Comput 67(2):791---794, 2016). In this paper, we have extended the stability analysis to the advection---diffusion equation, demonstrating that the tensor product FR formulation is stable on Cartesian meshes for the advection---diffusion case as well. The analysis in this paper also provides additional insights into the dependence on the VCJH parameter of the diffusion and stability characteristics of these schemes. Several numerical experiments that support the theoretical results are included.

A high-order semi-explicit discontinuous Galerkin solver for 3D incompressible flow with application to DNS and LES of turbulent channel flow

We present an efficient discontinuous Galerkin scheme for simulation of the incompressible Navier–Stokes equations including laminar and turbulent flow. We consider a semi-explicit high-order velocity-correction method for time integration as well as nodal equal-order discretizations for velocity and pressure. The non-linear convective term is treated explicitly while a linear system is solved for the pressure Poisson equation and the viscous term. The key feature of our solver is a consistent penalty term reducing the local divergence error in order to overcome recently reported instabilities in spatially under-resolved high-Reynolds-number flows as well as small time steps. This penalty method is similar to the grad–div stabilization widely used in continuous finite elements. We further review and compare our method to several other techniques recently proposed in literature to stabilize the method for such flow configurations. The solver is specifically designed for large-scale computations through matrix-free linear solvers including efficient preconditioning strategies and tensor-product elements, which have allowed us to scale this code up to 34.4 billion degrees of freedom and 147,456 CPU cores. We validate our code and demonstrate optimal convergence rates with laminar flows present in a vortex problem and flow past a cylinder and show applicability of our solver to direct numerical simulation as well as implicit large-eddy simulation of turbulent channel flow at Reτ=180 as well as 590.

Serendipity and Tensor Product Affine Pyramid Finite Elements

Using the language of finite element exterior calculus, we define two families of H 1 -conforming finite element spaces over pyramids with a parallelogram base. The first family has matching polynomial traces with tensor product elements on the base while the second has matching polynomial traces with serendipity elements on the base. The second family is new to the literature and provides a robust approach for linking between Lagrange elements on tetrahedra and serendipity elements on affinely-mapped cubes while preserving continuity and approximation properties. We define shape functions and degrees of freedom for each family and prove unisolvence and polynomial reproduction results.

A numerical analysis of the nodal Discontinuous Galerkin scheme via Flux Reconstruction for the advection-diffusion equation

Abstract A von Neumann analysis of the Flux Reconstruction (FR) formulation is performed for the linear advection-diffusion equation to investigate the stability, dissipation and dispersion associated with the nodal Discontinuous Galerkin (DG) scheme. We show that the maximum stable time step for advection-diffusion is stricter than that for pure-advection or pure-diffusion individually. However, the simple harmonic sum of the maximum stable time steps for advection and diffusion provides a suitable estimate for the linear advection-diffusion and Navier–Stokes equations on unstructured, tensor product elements. The estimates are accurate within 50% error on all test cases and are conservative on tests with Cartesian grids but not always on unstructured grids. We show that the CFL limit is strongly influenced by the choice of interface fluxes and, in general, the limit for a scheme using centered values is much higher than that which has one-sided values. We also verify that schemes with centered values produce less error for well resolved solutions while schemes with one-sided values produce less error for solutions that are under-resolved.

Hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes

We consider the hp-version interior penalty discontinuous Galerkin finite element method (DGFEM) for the numerical approximation of the advection-diffusion-reaction equation on general computational meshes consisting of polygonal/polyhedral (polytopic) elements. In particular, new hp-version a priori error bounds are derived based on a specific choice of the interior penalty parameter which allows for edge/face-degeneration. The proposed method employs elemental polynomial bases of total degree p (P_p-basis) defined in the physical coordinate system, without requiring the mapping from a given reference or canonical frame. Numerical experiments highlighting the performance of the proposed DGFEM are presented. In particular, we study the competitiveness of the p-version DGFEM employing a P_p-basis on both polytopic and tensor-product elements with a (standard) DGFEM employing a (mapped) Q_p-basis. Moreover, a computational example is also presented which demonstrates the performance of the proposed hp-version DGFEM on general agglomerated meshes.

Hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes

An hp-version interior penalty discontinuous Galerkin method (DGFEM) for the numerical solution of second-order elliptic partial differential equations on general computational meshes consisting of polygonal/polyhedral elements is presented and analyzed. Utilizing a bounding box concept, the method employs elemental polynomial bases of total degree p (𝒫p-basis) defined on the physical space, without the need to map from a given reference or canonical frame. This, together with a new specific choice of the interior penalty parameter which allows for face-degeneration, ensures that optimal a priori bounds may be established, for general meshes including polygonal elements with degenerating edges in two dimensions and polyhedral elements with degenerating faces and/or edges in three dimensions. Numerical experiments highlighting the performance of the proposed method are presented. Moreover, the competitiveness of the p-version DGFEM employing a 𝒫p-basis in comparison to the conforming p-version finite element method on tensor-product elements is studied numerically for a simple test problem.

An optimally blended finite-spectral element scheme with minimal dispersion for Maxwell equations

We study the dispersive properties of the time harmonic Maxwell equations for optimally blended finite-spectral element scheme using tensor product elements defined on rectangular grid in d-dimensions. We prove and give analytical expressions for the discrete dispersion relations for this scheme. We find that for a rectangular grid (a) the analytical expressions for the discrete dispersion error in higher dimensions can be obtained using one dimensional discrete dispersion error expressions; (b) the optimum value of the blending parameter is p/(p+1) for all p∈N and for any number of spatial dimensions; (c) analytical expressions for the discrete dispersion relations for finite element and spectral element schemes can be obtained when the value of blending parameter is chosen to be 0 and 1 respectively; (d) the optimally blended scheme guarantees two additional orders of accuracy compared with standard finite element and spectral element schemes; and (e) the absolute accuracy of the optimally blended scheme is O(p-2) and O(p-1) times better than that of the pure finite element and spectral element schemes respectively.

A spectral finite element approach to modeling soft solids excited with high-frequency harmonic loads

An approach for efficient and accurate finite element analysis of harmonically excited soft solids using high-order spectral finite elements is presented and evaluated. The Helmholtz-type equations used to model such systems suffer from additional numerical error known as pollution when excitation frequency becomes high relative to stiffness (i.e. high wave number), which is the case, for example, for soft tissues subject to ultrasound excitations. The use of high-order polynomial elements allows for a reduction in this pollution error, but requires additional consideration to counteract Runge’s phenomenon and/or poor linear system conditioning, which has led to the use of spectral element approaches. This work examines in detail the computational benefits and practical applicability of high-order spectral elements for such problems. The spectral elements examined are tensor product elements (i.e. quad or brick elements) of high-order Lagrangian polynomials with non-uniformly distributed Gauss–Lobatto–Legendre nodal points. A shear plane wave example is presented to show the dependence of the accuracy and computational expense of high-order elements on wave number. Then, a convergence study for a viscoelastic acoustic-structure interaction finite element model of an actual ultrasound driven vibroacoustic experiment is shown. The number of degrees of freedom required for a given accuracy level was found to consistently decrease with increasing element order. However, the computationally optimal element order was found to strongly depend on the wave number.