Truncation Error Estimation in the p-Anisotropic Discontinuous Galerkin Spectral Element Method

Abstract

In the context of discontinuous Galerkin spectral element methods (DGSEM), $\tau$-estimation has been successfully used for p-adaptation algorithms. This method estimates the truncation error of representations with different polynomial orders using the solution on a reference mesh of relatively high order. In this paper, we present a novel anisotropic truncation error estimator derived from the $\tau$-estimation procedure for the traditional DGSEM. We exploit the tensor product basis properties of the numerical solution to design a method where the total truncation error is calculated as a sum of its directional components. We show that the new error estimator is cheaper to evaluate than previous implementations of the $\tau$-estimation procedure and that it obtains more accurate extrapolations of the truncation error for representations of a higher order than the reference mesh. The robustness of the method allows performing the p-adaptation strategy with coarser reference solutions, thus further reducing the computational cost. The proposed estimator is validated using the method of manufactured solutions in a test case for the compressible Navier–Stokes equations.

Appendices

Appendix A: Isolated Truncation Error Dependence on Interpolation Error

According to Definition 5 and Eq. 37, the isolated truncation error in the DGSEM can be expressed for any basis function $\phi$ in an element e as

$$\begin{aligned} {\hat{\tau }}^N \bigr |_{{\varOmega }^e} = \hat{{\mathcal {R}}} ({\mathbf {I}}^N \bar{{\mathbf {q}}}) = \int ^N_{{\varOmega }^e} {\mathbf {s}}^N \phi d{\varOmega } + \int ^N_{{\varOmega }^e} {\mathscr {F}}^N \cdot \nabla \phi d {\varOmega } - \int ^N_{\partial {\varOmega }^e} {\mathscr {F}}^N \cdot {\mathbf {n}} \phi d \sigma , \end{aligned}$$

(71)

where the superindex N on the integrals indicates that they are approximated with a Gaussian quadrature of order N and the superindex e has been dropped for readability. Since the DGSEM is a collocation method, the value computed with Eq. 71 corresponds to the isolated truncation error on the node of the basis function $\phi$. The terms ${\mathbf {s}}^N$ and ${\mathscr {F}}^N$ can be expressed in terms of the interpolation error as

$$\begin{aligned} {\mathscr {F}}^N = {\mathbf {I}}^N {\mathscr {F}}(\bar{{\mathbf {q}}}) = {\mathscr {F}}(\bar{{\mathbf {q}}}) - \varepsilon ^N_{{\mathscr {F}}}, \ \ \ {\mathbf {s}}^N = {\mathbf {I}}^N {\mathbf {s}} = {\mathbf {s}} - \varepsilon ^N_{{\mathbf {s}}}. \end{aligned}$$

(72)

Inserting Eq. 72 into 71, integrating by parts, and expressing everything with $L_2({\varOmega })$ inner product notation we obtain,

(73)

where $(\cdot ,\cdot )^N_{{\varOmega }^e}$ stands for the $L_2$ product operator evaluated with a quadrature of order N in the domain ${{\varOmega }^e}$. The first term on the right-hand side vanishes since the value of $\varepsilon _{{\mathbf {s}}}^N$ is zero on the quadrature nodes (the DGSEM is a collocation method). Furthermore, it is reasonable to neglect the quadrature error since it is of a lower order of magnitude than the value of the integral. Therefore, we obtain

$$\begin{aligned} {\hat{\tau }}^N \bigr |_{{\varOmega }^e} \approx \left( \nabla \cdot \varepsilon ^N_{{\mathscr {F}}} , \phi \right) ^N_{{\varOmega }^e}. \end{aligned}$$

(74)

Appendix B: Anisotropic Non-isolated Truncation Error Estimation

In this section, we show briefly that the non-isolated truncation error can be estimated anisotropically using Proposition 1. In order to do so, we need some additional assumptions.

1.1 B.1 Additional Assumptions

As in Sect. 3.1, following assumptions are a consequence of the tensor product basis functions of the DGSEM and hold for sufficiently smooth solutions in the asymptotic range:

• (c) The discretization error has an anisotropic behavior and, therefore, can be decoupled in directional components. For the 2D case:

  $$\begin{aligned} \epsilon ^{N_1N_2} = \epsilon _1^{N_1N_2} + \epsilon _2^{N_1N_2}. \end{aligned}$$

  (75)

  As in (a), $\epsilon _i$ is the projection of the global discretization error, $\epsilon$, into a local direction, i.

• (d) The locally-generated discretization error in each direction depends only on the polynomial order in that direction:

  $$\begin{aligned} \epsilon _{{\varOmega },i}^{N_1N_2} = \epsilon _{{\varOmega },i}^{N_1N_2}(N_i) \end{aligned}$$

  (76)

Following the same reasoning as in Remark 2, and for reasons that will become clear at the end of the proof, an additional assumption is required:

• (e) The $\tau$-estimation procedure is performed element-wise while keeping the polynomial order in other elements sufficiently high so that:

  $$\begin{aligned} \left|\left|\epsilon ^N_{\partial {\varOmega }}\right|\right| \ll \left|\left|\epsilon ^N_{{\varOmega }}\right|\right| \end{aligned}$$

  (77)

As (a) and (b), assumptions (c) and (d) also follow from the work of Rubio et al. [42, 43]. The authors remark that assumptions (a), (b), (c) and (d) are consistent with the dependence of the non-isolated truncation error on the discretization error (Eq. 33). Furthermore, assumption (e) imposes an additional requirement for Proposition 1 to hold: that the polynomial orders of the elements not being analyzed must be maintained high enough during $\tau$-estimation, so that the externally-generated contributions to the discretization error are smaller than the internally-generated ones.

Let us note that the assumption (d) implies that, for smooth solutions in the asymptotic range, the discretization error in one direction does not change considerably when the polynomial order in another direction is changed:

$$\begin{aligned}&\epsilon _j^{N_iP_j} \approx \epsilon _j^{P_iP_j}, \nonumber \\&\epsilon _i^{N_iP_j} \ne \epsilon _i^{P_iP_j}, \end{aligned}$$

(78)

with $i \ne j$, and $1 \le i,j \le 2$.

Proof

Following the same procedure as in “Appendix A”, according to Definition 4 and Eq. 31, the non-isolated truncation error in the DGSEM can be expressed for any basis function $\phi$ in an element e as

$$\begin{aligned} \tau ^N \bigr |_{{\varOmega }^e} = {\mathcal {R}}^N ({\mathbf {I}}^N \bar{{\mathbf {q}}}) = \int ^N_{{\varOmega }^e} {\mathbf {s}}^N \phi d{\varOmega } + \int ^N_{{\varOmega }^e} {\mathscr {F}}({\mathbf {I}}^N\bar{{\mathbf {q}}}) \cdot \nabla \phi d {\varOmega } - \int ^N_{\partial {\varOmega }^e} {\mathscr {F}}^{*}({\mathbf {I}}^N \bar{{\mathbf {q}}},{\mathbf {I}}^N \bar{{\mathbf {q}}}^{\underline{\ }},{\mathbf {n}}) \phi d \sigma ,\nonumber \\ \end{aligned}$$

(79)

where $\bar{{\mathbf {q}}}^{\underline{\ }}$ is the external (neighbor element’s) solution and the superindex “e” has been dropped for the local solution. Since the DGSEM is a collocation method, the value computed with Eq. 79 corresponds to the non-isolated truncation error on the node of the basis function $\phi$. After inserting the definition of discretization error (Definition 2), $\bar{{\mathbf {q}}}=\bar{{\mathbf {q}}}^N+\epsilon ^N$, and expanding the fluxes using Taylor series we obtain

$$\begin{aligned} \tau ^N \bigr |_{{\varOmega }^e} \approx \int ^N_{{\varOmega }^e} \frac{\partial {\mathscr {F}}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N} \epsilon ^N \cdot \nabla \phi d {\varOmega } - \int ^N_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} \epsilon ^N \phi d \sigma - \int ^N_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}^{\underline{\ }}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} \epsilon ^{\underline{N}} \phi d \sigma ,\nonumber \\ \end{aligned}$$

(80)

where the interpolant of the discretization error is omitted for readability (${\mathbf {I}}^N\epsilon ^N \rightarrow \epsilon ^N$), $\epsilon ^N$ is the discretization error of the element e, and $\epsilon ^{\underline{N}}$ is the discretization error of a neighbor element connected through the surface $\partial {\varOmega }$. Notice that, for the sake of readability, the symbol for the external polynomial orders is the same as of the internal ones, i.e. N, although they can be different.

We now want to approximate the non-isolated truncation error through $\tau$-estimation. We part from the definition of the discretization error (Eq. 20). Adding and subtracting the discrete solution on a higher order grid, ${\mathbf {q}}^P$, yields

$$\begin{aligned} \epsilon ^N&= \bar{{\mathbf {q}}}-\bar{{\mathbf {q}}}^P+\bar{{\mathbf {q}}}^P-\bar{{\mathbf {q}}}^N \\ \epsilon ^N&= \epsilon ^P +\bar{ {\mathbf {q}}}^P -\bar{ {\mathbf {q}}}^N. \end{aligned}$$

Reorganizing we have

$$\begin{aligned} \bar{{\mathbf {q}}}^P = \bar{{\mathbf {q}}}^N + \epsilon ^N - \epsilon ^P. \end{aligned}$$

(81)

Therefore, the $\tau$-estimation yields

$$\begin{aligned} \tau _P^N \bigr |_{{\varOmega }^e}= & {} {\mathcal {R}} ({\mathbf {I}}^N \bar{{\mathbf {q}}}^P) \approx \int ^N_{{\varOmega }^e} \frac{\partial {\mathscr {F}}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N} (\epsilon ^N - \epsilon ^P) \cdot \nabla \phi d {\varOmega } - \int ^N_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} (\epsilon ^N - \epsilon ^P) \phi d \sigma \nonumber \\&- \int ^N_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}^{\underline{\ }}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} (\epsilon ^{\underline{N}} - \epsilon ^{\underline{P}}) \phi d \sigma . \end{aligned}$$

(82)

Since it is possible to decouple the discretization error inside our analyzed element in a locally-generated and an externally-generated component (Eq. 21), Eq. 82 can be rewritten as

$$\begin{aligned} \tau _P^N \bigr |_{{\varOmega }^e}&\approx \int ^N_{{\varOmega }^e} \frac{\partial {\mathscr {F}}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N} (\epsilon _{{\varOmega }}^N - \epsilon _{{\varOmega }}^P) \cdot \nabla \phi d {\varOmega } - \int ^N_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} (\epsilon _{{\varOmega }}^N - \epsilon _{{\varOmega }}^P) \phi d \sigma \nonumber \\&\quad + \int ^N_{{\varOmega }^e} \frac{\partial {\mathscr {F}}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N} (\epsilon _{\partial {\varOmega }}^N - \epsilon _{\partial {\varOmega }}^P) \cdot \nabla \phi d {\varOmega } - \int ^N_{\partial {\varOmega }^e}\frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} (\epsilon _{\partial {\varOmega }}^N - \epsilon _{\partial {\varOmega }}^P) \phi d \sigma \nonumber \\&\quad - \int ^N_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}^{\underline{\ }}} \biggr |_{\bar{{\mathbf {q}}}^N,\bar{{\mathbf {q}}}^{\underline{N}},{\mathbf {n}}} (\epsilon ^{\underline{N}} - \epsilon ^{\underline{P}}) \phi d \sigma . \end{aligned}$$

(83)

Equation 83 holds even for anisotropic representations, i.e. $N=(N_1,N_2,N_3)$ and $P=(P_1,P_2,P_3)$. Let us note that if the polynomial order of the elements that are not being analyzed is maintained as high as in the reference mesh, $\epsilon _{\partial {\varOmega }}^P$ cancels out $\epsilon _{\partial {\varOmega }}^N$ and $\epsilon ^{\underline{N}} - \epsilon ^{\underline{P}} \approx 0$, i.e., the $\tau$-estimation provides the locally-generated truncation error.

Let us now consider the case of 2D anisotropic coarsening in the direction i (${\mathbf {N}}=(N_i,P_j)$, ${\mathbf {P}}=(P_i,P_j)$). Taking into account assumptions (c) and (d), we obtain

$$\begin{aligned} \tau _{{\mathbf {P}}}^{{\mathbf {N}}} \bigr |_{{\varOmega }^e}&\approx \int ^{{\mathbf {N}}}_{{\varOmega }^e} \frac{\partial {\mathscr {F}}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^{{\mathbf {N}}}} (\epsilon _{{\varOmega },i}^{{\mathbf {N}}} - \epsilon _{{\varOmega },i}^{{\mathbf {P}}}) \cdot \nabla \phi d {\varOmega } - \int ^{{\mathbf {N}}}_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^{{\mathbf {N}}},\bar{{\mathbf {q}}}^{\underline{{{\mathbf {N}}}}},{\mathbf {n}}} (\epsilon _{{\varOmega },i}^{{\mathbf {N}}} - \epsilon _{{\varOmega },i}^{{\mathbf {P}}}) \phi d \sigma \nonumber \\&\quad + \int ^{{\mathbf {N}}}_{{\varOmega }^e} \frac{\partial {\mathscr {F}}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^{{\mathbf {N}}}} (\epsilon _{\partial {\varOmega }}^{{\mathbf {N}}} - \epsilon _{\partial {\varOmega }}^{{\mathbf {P}}}) \cdot \nabla \phi d {\varOmega } - \int ^{{\mathbf {N}}}_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}} \biggr |_{\bar{{\mathbf {q}}}^{{\mathbf {N}}},\bar{{\mathbf {q}}}^{\underline{{{\mathbf {N}}}}},{\mathbf {n}}} (\epsilon _{\partial {\varOmega }}^{{\mathbf {N}}} - \epsilon _{\partial {\varOmega }}^{{\mathbf {P}}}) \phi d \sigma \nonumber \\&\quad - \int ^{{\mathbf {N}}}_{\partial {\varOmega }^e} \frac{\partial \mathscr {F^*}}{\partial {\mathbf {q}}^{\underline{\ }}} \biggr |_{\bar{{\mathbf {q}}}^{{\mathbf {N}}},\bar{{\mathbf {q}}}^{\underline{{{\mathbf {N}}}}},{\mathbf {n}}} (\epsilon ^{\underline{{{\mathbf {N}}}}} - \epsilon ^{\underline{{{\mathbf {P}}}}}) \phi d \sigma . \end{aligned}$$

(84)

It is important to note that for sufficiently smooth solutions, the discretization errors in the high order mesh of order ${\mathbf {P}}$ are smaller than the discretization errors in the analyzed mesh ${\mathbf {N}}$ and therefore can be neglected. Finally, if assumptions (c) and (e) hold, the anisotropic version of Eq. 83 ($N=(N_1,N_2,N_3)$) can be reconstructed by summing all the directional components (Eq. 84) if the quadrature errors are neglected. I.e., we recover Eq. 45:

$$\begin{aligned} \tau ^{N_1N_2} \approx \tau ^{N_1P_2}_{P_1P_2} + \tau ^{P_1N_2}_{P_1P_2} \end{aligned}$$

(85)

$\square$

Appendix C: The Navier–Stokes Equations

The compressible Navier–Stokes equations in conservative form can be written in non-dimensional form as

$$\begin{aligned} {\mathbf {q}}_t + \nabla \cdot \left( {\mathscr {F}}^a - {\mathscr {F}}^{\nu } \right) = {\mathbf {s}}, \end{aligned}$$

(86)

where the conserved variables are ${\mathbf {q}} = \left( \rho , \rho u, \rho v, \rho w, \rho e \right) ^T$ ($\rho$ is the density; u, v and w are the velocity components; and e is the specific total energy), ${\mathbf {s}}$ is an external source term, and ${\mathscr {F}}^a$ and ${\mathscr {F}}^{\nu }$ are called the advective and diffusive flux dyadic tensors, respectively, which depend on ${\mathbf {q}}$. Expanding the fluxes in Cartesian coordinates leads to the expression,

$$\begin{aligned} {\mathbf {q}}_t + {\mathbf {f}}^a_x + {\mathbf {g}}^a_y + {\mathbf {h}}^a_z - \frac{1}{{\mathrm{Re}}} \left( {\mathbf {f}}^{\nu }_x + {\mathbf {g}}^{\nu }_y + {\mathbf {h}}^{\nu }_z \right) = {\mathbf {s}}. \end{aligned}$$

(87)

Here, Re is the Reynolds number. The advective fluxes are then defined as

$$\begin{aligned} {\mathbf {f}}^a = \begin{bmatrix} \rho u \\ p + \rho u^2\\ \rho u v \\ \rho u w \\ u (\rho e + p) \end{bmatrix}, {\mathbf {g}}^a = \begin{bmatrix} \rho v \\ \rho u v \\ p + \rho v^2\\ \rho v w \\ v (\rho e + p) \end{bmatrix}, {\mathbf {h}}^a = \begin{bmatrix} \rho w \\ \rho u w \\ \rho v w \\ p + \rho w^2\\ w (\rho e + p) \end{bmatrix}, \end{aligned}$$

(88)

where the pressure p is computed using the calorically perfect gas approximation. On the other hand, the diffusive fluxes are defined as

$$\begin{aligned} {\mathbf {f}}^{\nu }&= \begin{bmatrix} 0 \\ \tau _{xx} \\ \tau _{xy} \\ \tau _{xz} \\ u \tau _{xx} + v \tau _{xy} + w \tau _{xz} + \frac{\kappa }{(\gamma - 1) {\mathrm{Pr}} \mathrm {M}^2} T_x \end{bmatrix}, \end{aligned}$$

(89)

$$\begin{aligned} {\mathbf {g}}^{\nu }&= \begin{bmatrix} 0 \\ \tau _{yx} \\ \tau _{yy} \\ \tau _{yz} \\ u \tau _{yx} + v \tau _{yy} + w \tau _{yz} + \frac{\kappa }{(\gamma - 1) {\mathrm{Pr}} \mathrm {M}^2} T_y \\ \end{bmatrix}, \end{aligned}$$

(90)

$$\begin{aligned} {\mathbf {h}}^{\nu }&= \begin{bmatrix} 0 \\ \tau _{zx} \\ \tau _{zy} \\ \tau _{zz} \\ u \tau _{zx} + v \tau _{zy} + w \tau _{zz} + \frac{\kappa }{(\gamma - 1) {\mathrm{Pr}} \mathrm {M}^2} T_z \end{bmatrix}, \end{aligned}$$

(91)

where T is the temperature, $\gamma$ is the heat capacity ratio, and $\kappa$ is the thermal diffusivity. The nondimensional parameters are Pr, the Prandtl number; and M, the Mach number. The stress tensor components are computed using the Stokes hypothesis,

$$\begin{aligned} \tau _{ij}&= \mu \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) , i \ne j \end{aligned}$$

(92)

$$\begin{aligned} \tau _{ii}&= 2 \mu \left( \frac{\partial v_i}{\partial x_i} + \nabla \cdot {\mathbf {V}} \right) , \end{aligned}$$

(93)

with $\mu$ the fluid’s viscosity, and ${\mathbf {V}}$ the flow velocity. For the simulations in this paper we chose the typical parameters for air: ${\mathrm{Pr}} = 0.72$, $\gamma = 1.4$, while $\mu$ and $\kappa$ are calculated using Sutherland’s law.