Mimetic Metrics for the DGSEM

TL;DR

This paper introduces a new method to compute metric terms in DGSEMs that ensures they are divergence-free, which is important for accurate numerical simulations on curvilinear grids.

Contribution

A novel mimetic approach is introduced to compute divergence-free metric terms for DGSEMs using projections within de Rham Cohomology.

Findings

The new method guarantees divergence-free metric terms for DGSEMs.

The approach is based on projections that fit within the de Rham Cohomology framework.

Abstract

Free-stream preservation is an essential property for numerical solvers on curvilinear grids. Key to this property is that the metric terms of the curvilinear mapping satisfy discrete metric identities, i.e., have zero divergence. Divergence-free metric terms are furthermore essential for entropy stability on curvilinear grids. We present a new way to compute the metric terms for discontinuous Galerkin spectral element methods (DGSEMs) that guarantees they are divergence-free. The proposed mimetic approach uses projections that fit within the de Rham Cohomology.

Figures

Diagram of the continuous and discrete 3D de Rham complexFig. 2Diagram of the continuous and discrete 1D de Rham complex

Diagram of the continuous and discrete 2D de Rham complex

Free-stream preservation errors of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho e$$\end{document}$ with method of Kopriva [[4](#CR4)], Eq. (9), and the Mimetic approach, Eq. (19)Fig. 5Error in the metric terms with method of Kopriva, (9), and the Mimetic approach, see (19)

Full text

Introduction

We focus on discontinuous Galerkin spectral element methods (DGSEMs) for systems of conservation laws of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial }{\partial t} u_i + \text {div}_x\left( {\textbf {F}}_i(\textbf{u})\right) = 0 \qquad i = 1,\dots ,m \end{aligned}$$\end{document}

on a domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega \times [0,T]$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u} = \left( u_1, \dots , u_m\right) ^T$$\end{document} is the vector of conserved variables, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {F}}_i({\textbf {u}})$$\end{document} is the vector-valued flux function of the i-th equation and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {div}_x$$\end{document} denotes the divergence operator in physical space. In this paper, we denote differential operators in physical space with a subscript x, as in “div \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_x$$\end{document} ”. Furthermore, we use a bold font for vector-valued variables and functions, while we do not use bold letters for scalar quantities. In the following, we will only look at one of the equations at a time, so we will omit the subscript, i.

To solve the system numerically we partition the domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document} with a conforming mesh composed of quadrilateral (in 2D) or hexahedral (in 3D) elements.

Let

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} {\textbf {x}}: E & \longrightarrow Q, \\ \varvec{\xi } = \left( \xi ^1, \xi ^2, \xi ^3\right) ^T & \longmapsto {\textbf {x}} = \left( x_1, x_2, x_3\right) ^T \end{matrix}} \end{aligned}$$\end{document}

be the mapping from the reference element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E = [-1,1]^3$$\end{document} to an element Q in the physical domain.

For clarity, we will omit the subscript x when differential operators are applied with respect to the reference element coordinates.

After (1) is transformed to the reference element, it becomes

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial }{\partial t} \left( Ju\right) + \text {div}\left( (J {\textbf {a}}^s \cdot {\textbf {F}}(\textbf{u}))_{s=1}^3\right) = 0, \end{aligned}$$\end{document}

where J is the mapping Jacobian and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J{\textbf {a}}^i = J {\textbf {grad}}_x (\xi ^i)$$\end{document} are the metric terms,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} J{\textbf {a}}^i & = \frac{\partial {{\textbf {x}}}}{\partial {\xi ^j} } \times \frac{\partial {{\textbf {x}}}}{\partial {\xi ^k} } ,\,\, (i,j,k)\text { cyclic}, \\ \text { or equivalently [4]}, J{\textbf {a}}^i_n & = - \hat{{\textbf {x}}}_i \cdot {\textbf {curl}} \left( x_m {\textbf {grad}} \left( x_l\right) \right) ,\,\, (n,m,l)\text { cyclic}, \end{matrix}} \end{aligned}$$\end{document}

Here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\textbf {x}}}_i$$\end{document} are the canonical space unit vectors.

If we consider a constant state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}=\textbf{c}$$\end{document} and its corresponding flux \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {C}} = {\textbf {F}}(\textbf{c}) = const$$\end{document} , with constant or periodic boundary conditions, we get from (3)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial }{\partial t} \left( Ju\right) = - \text {div}\left( \left( J{\textbf {a}}^s \cdot {\textbf {C}}\right) _{s = 1}^3\right) = 0, \end{aligned}$$\end{document}

using that for the analytical global mapping the divergence of the metric terms is zero - the constant state stays constant in time.

For the DGSEM (see, e.g., [4]), we use a tensor product ansatz and interpolate Ju and the flux functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J {\textbf {a}}^s \cdot {\textbf {F}}(\textbf{u})$$\end{document} with polynomials of degree N in each direction. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I^N$$\end{document} be the tensor product interpolation operator and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l_{{\textbf {k}}}$$\end{document} the tensor product Lagrange polynomials. We multiply (3) with test functions, using the same Lagrange polynomials, and integrate over the reference element. Then we use integration-by-parts, replace the boundary flux with the numerical flux, use the summation-by-parts property and collocate interpolation and quadrature using Legendre-Gauss-Lobatto (LGL) nodes to arrive at the strong form of the DGSEM,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} & \frac{\partial }{\partial t} \left( Ju \right) _{{\textbf {i}}} + \text {div}\left( \left( I^N(J {\textbf {a}}^s \cdot {\textbf {F}}({\textbf {u}})) \right) _{s = 1}^3 \right) _{ {\textbf {i}}}\\ = & \frac{1}{\omega _{{\textbf {i}}}} \int _{\partial E,N} \left( I^N\left( F^{\text {num}}(\textbf{u}_l, \textbf{u}_r, {\textbf {n}})\right) - \left( I^N\left( {\textbf {F}} \cdot {\textbf {n}}\right) \right) \right) l_{{\textbf {i}}} \text { d}\xi , \end{matrix}} \end{aligned}$$\end{document}

for all multi-indices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {i}} \in \{0, \dots , N\}^n$$\end{document} . Here the multiindex \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {i}}$$\end{document} denotes the value of the function at the corresponding grid point, E is the reference element, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{{\textbf {i}}}$$\end{document} is the product of the LGL quadrature weights at the corresponding grid point and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _{\partial E, N}$$\end{document} denotes the numerical LGL-quadrature of the boundary integral (see [4]). Finally, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {n}}$$\end{document} denotes the outward directed normal vector on the interface and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_l$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_r$$\end{document} denote the left- and right-sided value on the interface, where left always means the value from ‘inside’ the element under consideration.

Substituting the constant state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{u}=\textbf{c}$$\end{document} into the DGSEM, the boundary term cancels due to the consistency of the numerical flux leaving

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial }{\partial t} \left( Ju \right) _{{\textbf {i}}} = - \text {div} \left( \left( I^N\left( J {\textbf {a}}^s \cdot {\textbf {C}}\right) \right) _{s=1}^3\right) _{ {\textbf {i}}}, \end{aligned}$$\end{document}

which is zero if the divergence of the discrete metric terms is zero,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {div} \left( \left( I^N\left( J{\textbf {a}}_i^s \right) \right) _{s = 1}^3 \right) = 0 \hspace{1cm} \forall i \in \{1,\dots ,3\}. \end{aligned}$$\end{document}

Even if the mapping \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {x}}$$\end{document} is polynomial, the metric terms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J{\textbf {a}}^i_n$$\end{document} computed using (4) are of a higher polynomial order. For the operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I^N$$\end{document} to still act as the identity operator in (8), we would have to restrict the mapping order to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \le \frac{N}{2}$$\end{document} . A different approach is used in [4] where the metric terms are computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} J{\textbf {a}}^i_n \approx -\frac{1}{2}\hat{{\textbf {x}}}_i \cdot {\textbf {curl}}\left( \textbf{I}^N\left( {x}_l {\textbf {grad}}\left( {x}_m\right) - {x}_m {\textbf {grad}}\left( {x}_l\right) \right) \right) \qquad (n,m,l) \text { cyclic}. \end{aligned}$$\end{document}

Here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{I}^N$$\end{document} is the component-wise interpolation. Since (9) takes the analytical curl of the interpolation and the divergence of a curl is zero, the curl form approximation of the metric terms guarantees free-stream preservation.

Mimetic Approach

We propose an alternative to the curl form of [4] using mimetic projections coming from finite element exterior calculus [1]. To apply this framework, we look at the de Rham cochain complex in 3D, and need to choose for each space a finite dimensional subspace \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_i$$\end{document} and a projection \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^i$$\end{document} so that the diagram 1 commutes.

To choose spaces and projections that fit our needs, we follow the work of Gerritsma et al. [2, 3] starting from the 1D de Rham complex 2.

We choose for each element the spaces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1 = \mathbb {P}^{N}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2 = \mathbb {P}^{N-1}$$\end{document} , which are the polynomial spaces of degree N and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N-1$$\end{document} respectively. To construct the projections we impose a sub-grid of the LGL nodes in each element, which gives an irregular distribution of sub-intervals between those nodes. Then we define basis functions on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} , where the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} basis corresponds to the LGL nodes and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} basis corresponds to the sub-intervals. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} these basis functions are just the Lagrange polynomials \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(l_i)_{i = 0}^N$$\end{document} , while for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} they are the so-called edge polynomials, e.g. [2], \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(h_i)_{i = 1}^N$$\end{document} defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} h_i(\xi ) = -\sum _{j = 0}^{i-1} \frac{\partial {l_j}}{\partial {\xi } } (\xi ), \end{aligned}$$\end{document}

with the property

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{\xi _{j-1}}^{\xi _j} h_i(\xi ) \text {d}\xi = \delta _{i,j}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{i,j}$$\end{document} is the Kronecker delta. Eq. (11) is a property analogous to that of the Lagrange polynomials, but for the sub-interval integrals instead of the point values. Using these basis functions we can define in each element the projections \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^2$$\end{document} by setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^1 = I^N$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^2$$\end{document} to be the histopolation, which is the approximation by polynomials preserving the sub-interval integrals. For one element this means

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p^2(f) = \sum _{i = 1}^N \left( \int _{\xi _{i-1}}^{\xi _i} f(\xi ) \text {d}\xi \right) h_i. \end{aligned}$$\end{document}

As shown in [2], the 1D diagram 2 commutes with these definitions. Since the LGL nodes include the boundary points of the element, we get inter-element continuity for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} while there is no such requirement for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} .Fig. 1. Diagram of the continuous and discrete 3D de Rham complexFig. 2Diagram of the continuous and discrete 1D de Rham complex

For the 3D case we can now use a tensor product ansatz of Lagrange and edge basis functions, and define the spaces as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V_1&= \mathbb {P}^{N} \otimes \mathbb {P}^{N} \otimes \mathbb {P}^{N},\end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V_2&= \left( \mathbb {P}^{N-1} \otimes \mathbb {P}^{N} \otimes \mathbb {P}^{N}, \mathbb {P}^{N} \otimes \mathbb {P}^{N-1} \otimes \mathbb {P}^{N}, \mathbb {P}^{N} \otimes \mathbb {P}^{N} \otimes \mathbb {P}^{N-1}\right) ^T,\end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V_3&= \left( \mathbb {P}^{N} \otimes \mathbb {P}^{N-1} \otimes \mathbb {P}^{N-1}, \mathbb {P}^{N-1} \otimes \mathbb {P}^{N} \otimes \mathbb {P}^{N-1}, \mathbb {P}^{N-1} \otimes \mathbb {P}^{N-1} \otimes \mathbb {P}^{N}\right) ^T,\end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V_4&= \mathbb {P}^{N-1} \otimes \mathbb {P}^{N-1} \otimes \mathbb {P}^{N-1}, \end{aligned}$$\end{document}

with continuity requirements derived analogously as in the 1D case, i.e. where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {P}^{N}$$\end{document} indicates directions with inter-element continuity, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {P}^{N-1}$$\end{document} indicates directions that are generally discontinuous. We note that these continuity requirements couple the nodal degrees of freedom similar to a continuous finite element approach.

The projection operators are similarly constructed based on the tensor product approach. We give the first component of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}^2: H({\textbf {curl}})\rightarrow V_2$$\end{document} as an example,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {p}}^2_1\left( {\textbf {f}}\right) = \sum _{i=1}^N \sum _{j = 0}^N \sum _{k = 0}^N \left( \int _{\xi _{i-1}^1}^{\xi _{i}^1} f_1(s,\xi _j^2,\xi ^3_k) \text { ds}\right) h_i(\xi ^1)l_j(\xi ^2)l_k(\xi ^3). \end{aligned}$$\end{document}

This construction is one of many possible approaches, for example the one used by Nédélec for cubic elements [6]. The Nédélec spaces are the same, but the Nédélec elements use higher moment degrees of freedom instead of sub-grid nodes, leading to different projections (see the discussion by Gerritsma et al. in [3]).

We can now define our approximative metric terms. Since the 3D diagram 1 commutes, we have two options with which to construct the discrete metric terms that result in the same approximation (up to machine precision).

Option red: Project \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {curl}}\left( x_m {\textbf {grad}}\left( x_l\right) \right) $$\end{document} directly via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}^3$$\end{document} to get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {p}}^3 \left( \left( J{\textbf {a}}^i_n\right) _{i = 1}^3\right) = {\textbf {p}}^3({\textbf {curl}}\left( x_m {\textbf {grad}} \left( x_l\right) \right) )\in V_3. \end{aligned}$$\end{document}

Option blue: Project \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_m {\textbf {grad}}\left( x_l\right) $$\end{document} via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}^2$$\end{document} to get \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}^2\left( x_m {\textbf {grad}} \left( x_l\right) \right) \in V_2$$\end{document} . We compute the metric terms by applying the curl, i.e.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {p}}^3 \left( \left( J{\textbf {a}}^i_n\right) _{i = 1}^3\right) = {\textbf {curl}}\left( {\textbf {p}}^2\left( x_m {\textbf {grad}}\left( x_l\right) \right) \right) \in V_3. \end{aligned}$$\end{document}

We then get free-stream preservation by construction:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {div}\left( {\textbf {p}}^3\left( {\textbf {curl}} \left( x_m {\textbf {grad}}\left( x_l\right) \right) \right) \right) = \text {div}\left( {\textbf {curl}}\left( {\textbf {p}}^2\left( x_m {\textbf {grad}}\left( x_l\right) \right) \right) \right) = 0 \end{aligned}$$\end{document}

**Remark 1 ** We note that for a fixed mesh, computing the metric terms is a pre-processing step with typically negligible impact on the overall compute time for 2D and 3D simulations. To make life simple, we could for instance first approximate the geometry with a piece-wise \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^0$$\end{document} polynomial ansatz by interpolation and use the resulting approximative mapping to compute the metric terms. By using standard numerical integration and differentiation procedures, the computational complexity is comparable to the approach presented by Kopriva [4].

It is also possible to directly use the (analytical) mapping as is. If possible, one could integrate and differentiate exactly. As an alternative to exact evaluation of metric terms, one could use again a numerical integration (and differentiation) approach that is implemented so that the errors are of machine precision magnitude, e.g., by using sufficient quadrature points and automatic differentiation.

**Remark 2 **In 2D, the approach in [4] is to interpolate the mapping components and then take the analytic 2D curl to achieve free-stream preservation. This approach is equivalent to the mimetic approach because the metric terms in 2D depend linearly on the mapping components, i.e., (4) reduces to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} J{\textbf {a}}^1 = \left( \frac{\partial {}}{\partial {\xi ^2} } x_2, - \frac{\partial {}}{\partial {\xi ^2} } x_1 \right) ^T,\qquad J{\textbf {a}}^2 = \left( - \frac{\partial {}}{\partial {\xi ^1} } x_2, \frac{\partial {}}{\partial {\xi ^1} } x_1 \right) ^T. \end{aligned}$$\end{document}

In other words, for 2D we have the smaller commuting diagram 3 where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {curl}}_v(f):= \left( \frac{\partial {f}}{\partial {\xi ^2} } , - \frac{\partial {f}}{\partial {\xi ^1} } \right) ^T$$\end{document} . The approach for 2D in [4] is equivalent to taking the blue path.Fig. 3. Diagram of the continuous and discrete 2D de Rham complex

Numerical Example

As an example, we discretize the 3D compressible Euler equations

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\ \rho e \end{pmatrix}_t + \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ \rho v_1 v_3 \\ \left( \rho e + p\right) v_1 \end{pmatrix}_{x_1} + \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ \rho v_2 v_3 \\ \left( \rho e + p\right) v_2 \end{pmatrix}_{x_2} + \begin{pmatrix} \rho v_3 \\ \rho v_1 v_3 \\ \rho v_2 v_3 \\ \rho v_3^2 + p\\ \left( \rho e + p\right) v_3 \end{pmatrix}_{x_3} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \end{aligned}$$\end{document}

where e is the specific total energy and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p = \left( \gamma - 1\right) \left( \rho e - \frac{1}{2} \rho \left( v_1^2 + v_2^2 + v_3^2 \right) \right) , \end{aligned}$$\end{document}

with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 1.4$$\end{document} and initial condition

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho \equiv 1.0\text {, }{ \rho v_1 \equiv 0.1} \text {, } \rho v_2 \equiv -0.2 \text {, } \rho v_3 \equiv 0.7 \text {, } \rho e \equiv 10.0. \end{aligned}$$\end{document}

We use the mapping

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {x}}(\varvec{\xi }) = \varvec{\xi } + \theta (\varvec{\xi }) (1,1,1)^T \end{aligned}$$\end{document}

from the reference domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-1,1]^3$$\end{document} with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta (\varvec{\xi }) := 0.1 \cos (\pi \xi ^1) \cos (\pi \xi ^2) \cos (\pi \xi ^3), \end{aligned}$$\end{document}

and end time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = 1$$\end{document} . We use a conforming mesh with two elements in each direction and periodic boundary conditions. For time integration we use a fourth-order, five-stage low storage Runge Kutta method by Carpenter and Kennedy ([5], page 13, solution 3) with a CFL of 0.2, which is stable for all polynomial degrees investigated. The time step \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta t$$\end{document} is calculated by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varDelta t = \text {CFL} \frac{2}{|\lambda _{\text {max}} (N+1)|}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{\text {max}}$$\end{document} is the maximum wave speed of the PDE system evaluated at the LGL nodes scaled by the geometry. Simulations are done with a slightly modified version of the DGSEM code of the Julia package Trixi.jl (see [7] and [8]) performed on a Macbook Pro M2, single thread, with MacOS 14.5.

We examine both free-stream preservation and the approximation of the analytic metric terms by varying the polynomial degree up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 25$$\end{document} for both the curl form by Kopriva in [4] and the mimetic variant. Both errors are measured in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\infty }$$\end{document} norm. The norms are computed numerically using LGL integration with 51 points in each direction in every element for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} norm, and by the maximum of the point value errors on the same points as the discrete evaluation of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\infty }$$\end{document} norm. For the metric terms we compute on each of the 51 LGL points in each direction per element the norms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left| \left| \text {vec}\left( Ja_n^i\right) _{i,n = 1}^3\right| \right| _2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left| \left| \text {vec}\left( Ja_n^i\right) _{i,n = 1}^3\right| \right| _{\infty }$$\end{document} and then apply the LGL integration formula for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} or choose the maximum value for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\infty }$$\end{document} norm respectively.Fig. 4. Free-stream preservation errors of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho e$$\end{document} with method of Kopriva [4], Eq. (9), and the Mimetic approach, Eq. (19)Fig. 5. Error in the metric terms with method of Kopriva, (9), and the Mimetic approach, see (19)

In Fig. 4a and 4b we can see that free-stream preservation is achieved with both the curl and the mimetic approach, and that the rounding errors grow with increasing polynomial degree for both. We also observe that the mimetic approach has rounding errors smaller than Kopriva’s curl form by over an order of magnitude in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\infty }$$\end{document} norm.

In Figure 5a and 5b we see convergence of the metric terms to machine precision, though the rounding errors are again higher for Kopriva’s approximation.

Conclusion

We have constructed an alternative way to compute the metric terms for the DGSEM scheme in a divergence-free way via a mimetic approach using de Rham complexes for both two and three space dimensions. In two space dimensions the method of [4] and the mimetic approach are analytically the same, but they differ in three dimensions. Both approaches do obtain free-stream preservation for conforming meshes and exhibit a similar convergence behavior for the metric terms. However, in all cases the mimetic approach has better rounding error properties. We remark that the finite element exterior calculus projections do not only exist for quadrilateral or hexahedral meshes but also, for example, for triangular or tetrahedral meshes [1].