Abstract
The flux reconstruction (FR) methodology provides a unifying description of many high-order schemes, including a particular discontinuous Galerkin (DG) scheme and several spectral difference (SD) schemes. In addition, the FR methodology has been used to generate new classes of high-order schemes, including the recently discovered ‘energy stable’ FR schemes. These schemes, which are often referred to as VCJH (Vincent–Castonguay–Jameson–Huynh) schemes, are provably stable for linear advection–diffusion problems in 1D and on triangular elements. The VCJH schemes have been successfully applied to a wide variety of problems in 1D and 2D, ranging from linear advection–diffusion problems, to fluid mechanics problems requiring the solution of the compressible Navier–Stokes equations. Based on the results of these numerical experiments, it has been shown that certain VCJH schemes maintain the expected order of spatial accuracy and possess explicit time-step limits which rival those of the collocation-based nodal DG scheme. However, it remained to be seen whether the VCJH schemes could be extended to 3D on tetrahedral elements, enabling their convenient application to the complex geometries that arise in many real-world problems. For the first time, this article presents an extension of the VCJH schemes to tetrahedral elements. This work provides a formal proof of the stability of the new schemes and assesses their performance via numerical experiments on model problems.
Similar content being viewed by others
References
Allaneau, Y., Jameson, A.: Connections between the filtered discontinuous Galerkin method and the flux reconstruction approach to high order discretizations. Comput. Methods Appl. Mech. Eng. 200(49), 3628–3636 (2011)
Arnold, D.N.: An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal. 19(4), 742–760 (1982)
Bassi, F., Rebay, S.: Accurate 2D Euler computations by means of a high order discontinuous finite element method. In: 14th International Conference on Numerical Methods in Fluid Dynamics. Bangalor, India (1994)
Bassi, F., Rebay, S.: Discontinuous finite element high order accurate numerical solution of the compressible Navier–Stokes equations. In: ICFD Conference on Numerical Methods in Fluid Dynamics. University of Oxford, Oxford, England (1995)
Bassi, F., Rebay, S.: A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations. J. Comput. Phys. 131(2), 267–279 (1997)
Bassi, F., Rebay, S., Mariotti, G., Pedinotti, S., Savini, M.: A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows. In: Decuypere, R., Dibelius, G. (eds.) 2nd European Conference on Turbomachinery Fluid Dynamics and Thermodynamics. Antwerpen, Belgium (1997)
Brady, M., Horn, B.K.P.: Rotationally symmetric operators for surface interpolation. Comput. Vis. Graph. Image Process. 22(1), 70–94 (1983)
Carpenter, M.H., Kennedy, C.: Fourth-Order 2N-Storage Runge–Kutta Schemes. Tech. Rep. TM 109112, NASA, Langley Research Center (1994)
Castonguay, P.: High-Order Energy Stable Flux Reconstruction Schemes for Fluid Flow Simulations on Unstructured Grids. Ph.D. Thesis, Stanford University (2012)
Castonguay, P., Vincent, P.E., Jameson, A.: A new class of high-order energy stable flux reconstruction schemes for conservation laws on triangular grids. J. Sci. Comput. 51(1), 224–256 (2011)
Cockburn, B., Hou, S., Shu, C.W.: The Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case. Math. Comput. 54(190), 545–581 (1990)
Cockburn, B., Shu, C.W.: The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal. 35(6), 2440–2463 (1998)
Danielsson, P.E., Seger, O.: Rotation invariance in gradient and higher order derivative detectors. Comput. Vis. Graph. Image Process. 49(2), 198–221 (1990)
Friedrichs, K.O.: Symmetric hyperbolic linear differential equations. Commun. Pure Appl. Math. 7, 345–392 (1954)
Gao, H., Wang, Z.J.: A high-order lifting collocation penalty formulation for the Navier–Stokes equations on 2D mixed grids. In: 19th AIAA Computational Fluid Dynamics. San Antonio, TX (2009)
Grimson, W.E.L.: From Images to Surfaces: A Computational Study of the Human Early Visual System. MIT Press, Cambridge (1981)
Haga, T., Gao, H., Wang, Z.J.: A high-order unifying discontinuous formulation for 3D mixed grids. In: 48th AIAA Aerospace Sciences Meeting. Orlando, FL (2010)
Haga, T., Gao, H., Wang, Z.J.: A high-order unifying discontinuous formulation for the Navier–Stokes equations on 3D mixed grids. Mathe. Model. Nat. Phenom. 6(3), 28–56 (2011)
Hesthaven, J.S., Warburton, T.: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer, Berlin (2007)
Hirsch, C.: Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics, Volume I. Butterworth-Heinemann, London (2007)
Hu, F.Q., Hussaini, M.Y., Rasetarinera, P.: An analysis of the discontinuous Galerkin method for wave propagation problems. J. Comput. Phys. 151(2), 921–946 (1999)
Huynh, H.T.: A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. In: 18th AIAA Computational Fluid Dynamics Conference. Miami, FL (2007)
Huynh, H.T.: A reconstruction approach to high-order schemes including discontinuous Galerkin for diffusion. In: 47th AIAA Aerospace Sciences Meeting. Orlando, FL (2009)
Huynh, H.T.: High-order methods including discontinuous Galerkin by reconstructions on triangular meshes. In: 49th AIAA Aerospace Sciences Meeting. Orlando, FL (2011)
Jameson, A.: A proof of the stability of the spectral difference method for all orders of accuracy. J. Sci. Comput. 45(1), 348–358 (2010)
Kannan, R., Wang, Z.J.: LDG2: a variant of the LDG flux formulation for the spectral volume method. J. Sci. Comput. 46(2), 314–328 (2011)
Kopriva, D.A., Kolias, J.H.: A conservative staggered-grid Chebyshev multidomain method for compressible flows. J. Comput. Phys. 125, 244–261 (1996)
Lax, P.D.: Weak solutions of nonlinear hyperbolic equations and their numerical computation. Commun. Pure Appl. Math. 7, 159–193 (1954)
Liu, Y., Vinokur, M., Wang, Z.J.: Spectral difference method for unstructured grids I: basic formulation. J. Comput. Phys. 216, 780–801 (2006)
Peraire, J., Persson, P.O.: The compact discontinuous Galerkin (CDG) method for elliptic problems. SIAM J. Sci. Comput. 30(4), 1806–1824 (2009)
Raviart, P., Thomas, J.: A mixed finite element method for second-order elliptic problems. In: Mathematical Aspects of Finite Element Methods, pp. 292–315. Springer, Berlin (1977)
Van den Abeele, K., Lacor, C.: An accuracy and stability study of the 2D spectral volume method. J. Comput. Phys. 226(1), 1007–1026 (2007)
Vincent, P.E., Castonguay, P., Jameson, A.: A new class of high-order energy stable flux reconstruction schemes. J. Sci. Comput. 47(1), 50–72 (2011)
Vincent, P.E., Jameson, A.: Facilitating the adoption of unstructured high-order methods amongst a wider community of fluid dynamicists. Math. Model. Nat. Phenom. 6(3), 97–140 (2011)
Wang, Z.J., Gao, H.: A unifying lifting collocation penalty formulation for the Euler equations on mixed grids. In: AIAA P. 47th AIAA Aerospace Sciences Meeting, Orlando, FL, Jan. 5–8 (2009)
Wang, Z.J., Gao, H.: A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. J. Comput. Phys. 228(21), 8161–8186 (2009)
Williams, D.M.: Energy Stable High-Order Methods for Simulating Unsteady, Viscous, Compressible Flows on Unstructured Grids. Ph.D. Thesis, Stanford University (2013)
Williams, D.M., Castonguay, P., Vincent, P.E., Jameson, A.: Energy stable flux reconstruction schemes for advection–diffusion problems on triangles. J. Comput. Phys. 250, 53–76 (2013)
Yu, M., Wang, Z.J.: On the connection between the correction and weighting functions in the correction procedure via reconstruction method. J. Sci. Comput. 54(1), 227–244 (2012)
Acknowledgments
The authors would like to thank the National Science Foundation Graduate Research Fellowship Program, the Stanford Graduate Fellowships program, the National Science Foundation (Grants 0708071 and 0915006), the Air Force Office of Scientific Research (Grants FA9550-07-1-0195 and FA9550-10-1-0418) and NVIDIA for supporting this work.
Author information
Authors and Affiliations
Corresponding author
Appendices
Appendix A: Constructing the Energy Stable (VCJH) Correction Fields
In this section, a procedure will be presented for constructing energy stable (VCJH) correction fields \(\phi _{f,l}\) and \(\psi _{f,l}\) that satisfy Eqs. (43) and (53), respectively.
1.1 Preliminaries
Recall that \(\phi _{f,l} \equiv \hat{\nabla } \cdot \mathbf{h}_{f,l}\) and \(\psi _{f,l} \equiv \hat{\nabla } \cdot \mathbf{g}_{f,l}\). In accordance with these definitions, it may appear natural to first construct precise expressions for the correction functions \(\mathbf{h}_{f,l}\) and \(\mathbf{g}_{f,l}\), and thereafter apply the divergence operator \(\hat{\nabla }\) to these expressions in order to obtain \(\phi _{f,l}\) and \(\psi _{f,l}\). However, this is not the best approach, as \(\mathbf{h}_{f,l}\) and \(\mathbf{g}_{f,l}\) may not be unique. In particular, for \(p > 1\) there are an unlimited number of correction functions which have the same divergence (or effectively, the same correction field). In addition, the implementation of the VCJH approach only requires the precise formulation of the correction fields \(\phi _{f,l}\) and \(\psi _{f,l}\), and not of the correction functions themselves. Specifically, the VCJH approach only requires definitions of the normal components of the correction functions (\(\mathbf{h}_{f,l} \cdot \hat{\mathbf{n}}\) and \(\mathbf{g}_{f,l} \cdot \hat{\mathbf{n}}\)) which are given by Eqs. (23) and (24). In what follows, it will be shown that Eqs. (23), (24), (43), and (53) are sufficient for constructing precise definitions of \(\phi _{f,l}\) and \(\psi _{f,l}\).
1.2 Constructing the Correction Fields \(\phi _{f,l}\)
One may arrive at a procedure for constructing the correction fields \(\phi _{f,l}\) by manipulating Eq. (43) and utilizing the definition of \(\mathbf{h}_{f,l} \cdot \hat{\mathbf{n}}\) (Eqs. (23) and (24)). Recall from section 3, Lemma 3.1, that the correction functions \(\mathbf{h}_{f,l}\) and correction fields \(\phi _{f,l}\) are required to satisfy Eq. (43), which can be expanded as follows
Upon rearranging and simplifying Eq. (144), one obtains the following
Here, each basis function \(\hat{\ell }_{j}\) is equal to the product of \(J_k\) and \(\ell _{j}\), each normal component of a correction function \((\mathbf{h}_{f,l} \cdot \hat{\mathbf{n}})\) is defined by Eq. (23) (with \(\mathbf{h}_{f,l}\) in place of \(\mathbf{g}_{f,l}\)), and each correction field \(\phi _{f,l}\) remains to be determined. In order to begin the process of computing \(\phi _{f,l}\), note that it is a degree \(p\) polynomial function which can be expressed as follows
where each \(\sigma _{\jmath }\) is a constant coefficient (yet to be computed) and each basis function \(L_{\jmath } \left( \hat{\mathbf{x}} \right) \) is a member of an orthonormal basis of degree \(p\) on the reference element \(\varvec{\Omega }_S\) defined as follows
where \(P_n^{\left( \alpha , \beta \right) }\) is the normalized \(n^{th}\) order Jacobi polynomial (as defined in [19]), and where
In addition, note that each function \(\hat{\ell }_{j}\) in Eq. (145) is a degree \(p\) polynomial which can be expressed as follows
where each constant coefficient \(\gamma _{\imath }\) is a known quantity because each function \(\hat{\ell }_j\) is a known quantity.
On substituting Eqs. (146) and (149) into Eq. (145), one obtains
or equivalently
Upon noting that \(L_{\imath }\) and \(L_{\jmath }\) are orthonormal polynomials, one may derive the following expression from Eq. (151)
Equation (152) holds for \(\imath = 1, \ldots , N_p\) and provides a system of \(N_p\) equations for the \(N_p\) unknown coefficients \(\sigma _{\jmath }\). Together, Eqs. (152) and (146) completely define \(\phi _{f,l}\).
1.3 Constructing the Correction Fields \(\psi _{f,l}\)
In following the approach of the previous section, one may arrive at a procedure for constructing the correction fields \(\psi _{f,l}\) by manipulating Eq. (53) and utilizing the definition of \(\mathbf{g}_{f,l} \cdot \hat{\mathbf{n}}\) (Eqs. (23) and (24)). Once these manipulations (which are omitted for the sake of brevity) are performed, one obtains the following formula for each field \(\psi _{f,l}\)
where the \(N_p\) unknown coefficients \(\xi _{\jmath }\) can be obtained from the following system of \(\imath = 1, \ldots , N_p\) equations
Appendix B: Energy Stable (VCJH) Filter Matrices
This section discusses the procedure for forming the energy stable (VCJH) filter matrices and examines the resulting structure of the filter matrices.
1.1 Procedure for Forming the Filter Matrices
The filters \(\mathbf{F}_1\) and \(\mathbf{F}_2\) are obtained from the matrices \(\mathbf{M}^{k}, \mathcal M ^{k}, \mathbf{K}^{k}\), and \(\mathcal K ^{k}\) via Eq. (142). The mass matrices \(\mathbf{M}^{k}\) and \(\mathcal M ^{k}\) can be constructed from inner products of the nodal basis functions \(\ell _j \left( \hat{\mathbf{x}} \right) \) in accordance with Eq. (70). However, it is often difficult to compute the nodal basis functions \(\ell _j \left( \hat{\mathbf{x}} \right) \) directly, and one is often better served by using the orthonormal basis functions \(L_{j} \left( \hat{\mathbf{x}} \right) \) in their place. In particular, one may use the orthonormal basis functions to define a ‘Vandermonde matrix’ \(\mathbf{V}^k\) which has the following entries
and next, one may use \(\mathbf{V}^k\) to compute the mass matrices as follows
Note that the derivation of Eqs. (156) and (157) is discussed in detail in [19].
In a similar fashion, one may construct expressions for \(\mathbf{K}^{k}\) and \(\mathcal K ^{k}\) in terms of the orthonormal basis functions. Recall that Eqs. (76) and (77) provide expressions for \(\mathbf{K}^{k}\) and \(\mathcal K ^{k}\) in terms of the matrices \(\left( \hat{\mathbf{D}}^{(p,v,w)} \right) ^{T} \hat{\mathbf{D}}^{(p,v,w)}\), where each \(\hat{\mathbf{D}}^{(p,v,w)}\) is a matrix which has the following entries
One may use the orthonormal basis functions \(L_{j} \left( \hat{\mathbf{x}} \right) \) to construct similar matrices \(\hat{\hat{\mathbf{D}}}^{(p,v,w)}\) which have the following entries
One may then relate each \(\hat{\hat{\mathbf{D}}}^{(p,v,w)}\) to each \(\hat{\mathbf{D}}^{(p,v,w)}\) by using the following identity
Note that the derivation of Eq. (160) is discussed in detail in [19].
Upon substituting Eq. (160) into Eqs. (76) and (77), one obtains the following expressions for \(\mathbf{K}^{k}\) and \(\mathcal K ^{k}\)
Next, on substituting Eqs. (161), (162), (156), and (157) into Eq. (142), one obtains the following expressions for \(\mathbf{F}_1\) and \(\mathbf{F}_2\)
where
and
where
and
Note that Eqs. (164) and (166) define new filtering matrices \(\hat{\hat{\mathbf{F}}}_1\) and \(\hat{\hat{\mathbf{F}}}_2\). These matrices can be viewed as filters which act on the orthonormal basis (in contrast to \(\mathbf{F}_1\) and \(\mathbf{F}_2\) which can be viewed as filters which act on the nodal basis). Conveniently, the two sets of filters are related via left and right multiplication by the Vandermonde matrix and its inverse (as shown in Eqs. (163) and (165)).
Now, having established a method for constructing \(\mathbf{F}_1\) and \(\mathbf{F}_2\) using the orthornormal basis (via \(\hat{\hat{\mathbf{F}}}_1\) and \(\hat{\hat{\mathbf{F}}}_2\)), one may obtain insights into the overall effects of the filtering process by examining the sparsity patterns of \(\hat{\hat{\mathbf{F}}}_1\) and \(\hat{\hat{\mathbf{F}}}_2\).
1.2 Sparsity Patterns of the Filter Matrices
On evaluating Eq. (164), one finds that \(\hat{\hat{\mathbf{F}}}_1\) has the following block structure
where \(\mathbf{I}_{1}^{B} \in \mathbb R ^{N_p^{\ell } \times N_p^{\ell }}\) is an identity matrix, \(\mathbf{F}_{1}^{B} \in \mathbb R ^{N_p^{u} \times N_p^{u}}\) is a dense matrix of filtering coefficients, \(N_p^{\ell } = N_p - N_p^{u}\) is the number of orthonormal basis functions of degree \(\le \left( p - 1\right) \), and \(N_p^{u} = \frac{1}{2} \left( p + 1\right) \left( p + 2 \right) \) is the number of orthonormal basis functions of degree \(p\). The structure of \(\hat{\hat{\mathbf{F}_1}}\) ensures that only the degree \(p\) orthonormal basis functions are effected by the filtering matrix. All basis functions of degree \(\le \left( p-1\right) \) are multiplied by the identity matrix and remain unaffected.
Similarly, on evaluating Eq. (166), one finds that \(\hat{\hat{\mathbf{F}}}_2\) has the following structure
where \(\mathbf{I}_{2}^{B} \in \mathbb R ^{N_p^{\ell } \times N_p^{\ell }}\) is an identity matrix and \(\mathbf{F}_{2}^{B} \in \mathbb R ^{N_p^{u} \times N_p^{u}}\) is a dense matrix of filtering coefficients. The structure of \(\hat{\hat{\mathbf{F}_2}}\) is similar to the structure of \(\hat{\hat{\mathbf{F}_1}}\), as the structure of \(\hat{\hat{\mathbf{F}}}_2\) also ensures that only the degree \(p\) orthonormal basis functions are effected by the filtering matrix.
In summary, the filters \(\mathbf{F}_1\) and \(\mathbf{F}_2\) are related (via the Vandermonde matrix) to the filters \(\hat{\hat{\mathbf{F}_1}}\) and \(\hat{\hat{\mathbf{F}_2}}\) which act on the orthonormal basis functions \(L_j (\hat{\mathbf{x}})\), and effect only the highest (degree \(p\)) modes of the residual and the auxiliary variable, respectively.
Rights and permissions
About this article
Cite this article
Williams, D.M., Jameson, A. Energy Stable Flux Reconstruction Schemes for Advection–Diffusion Problems on Tetrahedra. J Sci Comput 59, 721–759 (2014). https://doi.org/10.1007/s10915-013-9780-2
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-013-9780-2