Truncation Error Estimation in the pAnisotropic Discontinuous Galerkin Spectral Element Method
Abstract
In the context of Discontinuous Galerkin Spectral Element Methods (DGSEM), estimation has been successfully used for padaptation 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 estimation procedure for 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 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 padaptation 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 NavierStokes equations.
Keywords:
Highorder discontinuous Galerkin Spectral methods pAnisotropic representations Estimation Truncation error Anisotropic padaptationMsc:
65M15 65M50 65M60 65M70∎
1 Introduction
Highorder Discontinuous Galerkin (DG) methods are becoming a popular alternative to low order methods for solving Partial Differential Equations (PDEs) because of their high accuracy and flexibility Wang2013High (); Cockburn2000 (). Among those, the Discontinuous Galerkin Spectral Element Method (DGSEM) kopriva2009implementing (); Kopriva2002 () is a nodal (collocation) version of the DG method on hexahedral meshes which allows panisotropic representations and has been used in a wide range of applications Kopriva2002 (); Rasetarinera2001 (); Acosta2011 (); Deng2007 (); Fraysse2016 (). In the DG approach, the continuity constraint on element interfaces is relaxed, allowing for discontinuities in the numerical solution. This feature makes them more robust than continuous methods for describing advectiondominated problems, like the ones usually encountered in fluid dynamics. Moreover, DG methods can handle nonconforming meshes with hanging nodes and/or different polynomial orders efficiently, as is necessary for mesh adaptation strategies Riviere2008 (); Kopriva2002 (); Ferrer2012 ().
Error estimates are a powerful tool in computational sciences as they quantify how accurately a numerical solution satisfies the governing mathematical equations oberkampf2010verification (); roache1998verification (); Roy2010 (). A precise assessment of the numerical errors is useful for defect correction (a technique that enables high accuracy by correcting the numerical solution using an estimation of the error roache1998verification (); Martin1996 ()), or for guiding mesh adaptation strategies Zienkiewicz2006 (); Mavriplis1989 (); Mavriplis1994 (). The former requires highly accurate estimates of the discretization error and, therefore, a significant amount of computational resources is usually invested in computing them Phillips2011 (). The latter has been broadly studied in the literature. In particular, the most common approaches are the adjointbased adaptation Estep1995 (); Hartmann2002 (); Hartmann2006 (); Venditti2002 (); Pierce2004 (), where the numerical error of a functional (e.g. lift or drag) is estimated, which involves a high computational cost; the featurebased adaptation, which relies on on easytocompute adaptation criteria, such as the assessment of jumps across element interfaces in the case of DG discretizations Remacle2003 (), or the identication of large gradients Aftosmis1994 (); Persson2006 (); and the localerrorbased adaptation Mavriplis1989 (); Mavriplis1994 (); berger1987adaptive (); Fraysse2012 (); Syrakos2012 (); Syrakos2006 (); Kompenhans2016a (); Kompenhans2016 (), which depends on the assessment of any measurable local error in all the cells of the domain. A detailed comparison of the different approaches for error estimation and adaptation can be found in Fraysse2012 () in the context of finite volumes or Kompenhans2016a () for highorder DG schemes. The localerrorbased adaptation methods are interesting since, in contrast to featurebased methods, they provide a way to predict and control the overall accuracy, and are computationally cheaper than adjointbased schemes Kompenhans2016a (); Kompenhans2016 (). The topic of our work is the development of an accurate and cheap local error estimator to drive panisotropic adaptation in the DGSEM.
Two different errors are particularly relevant. On the one hand, the discretization error is the most important, but also the most difficult error to estimate Phillips2011 (). It is defined as the difference between the exact and numerical solutions to the problem and can be approximated by means of solving the Discretization Error Transport Equation (DETE) Acosta2011 (), an auxiliary PDE whose approximation involves the investment of further computational resources. Some of the first works using estimations of the local discretization error in highorder methods were proposed by Mavriplis Mavriplis1989 (); Mavriplis1994 (), who developed hpadaptation techniques for the Spectral Element Method, and Casoni et al. EvaCasoniyAntonioHuerta2011 (), who used a similar approach to evaluate where to add artificial viscosity for shock capturing in Discontinuous Galerkin discretizations.
On the other hand, the truncation error is defined as the difference between the discrete partial differential operator and the exact partial differential operator,
(1) 
and is usually evaluated for the exact solution of the PDE Phillips2011 (); Rubio2015 (); Kompenhans2016 (); Kompenhans2016a (). The truncation error is related to the discretization error through the DETE Roy2010 (), where it acts as a local source term. This relation makes it useful as an indicator for mesh adaptation methods Syrakos2012 (); Choudhary2013 () since refining the mesh where the truncation error is high reduces the discretization error in all the mesh Rubio2015 (), with an additional advantage: the truncation error estimation requires less computational effort. Furthermore, in hyperbolic problems the discretization error is strongly advected, i.e. it is transmitted downstream from underresolution areas, but the truncation error is only weakly advected. Therefore, an adaptation procedure based on the truncation error targets specifically the underresolved areas, whereas one based on the discretization error targets the underresolved areas and the zones downstream of them rubio2015truncation (); Rubio2015 (). This makes the truncation error more suitable for adaptation purposes than the discretization error. Finally, it has been shown that controlling the truncation error targets the numerical accuracy of all functionals at once Kompenhans2016 (), ensuring that adapting a mesh using the truncation error leads necessarily to an error decrease in any other functional (e.g. lift or drag). For all these reasons, we focus on truncation error estimators in this paper.
From a practical point of view, the truncation error can be estimated using a hierarchy of meshes. On the one hand, Venditti and Darmofal Venditti2002 () and Phillips et al. Phillips2017 (); Phillips2013 () studied the possibility of estimating the truncation error by evaluating a coarse grid solution in the partial differential operator of a fine grid, an approach known as the coarsetofine approach. On the other hand, the finetocoarse approach, also known as the estimation method, was introduced by Brandt Brandt1984 () and consists in estimating the local truncation error by using a fine grid solution interpolated to the coarse grid. Phillips Phillips2014 () showed that the finetocoarse (estimation) method produces more accurate results than the coarsetofine approach and, therefore, it is the one retained in this work.
The estimation approach has been successfully used for adaptation purposes in loworder Finite Difference berger1987adaptive () and Finite Volume schemes Fraysse2012 (); Syrakos2012 (); Syrakos2006 (). Moreover, Rubio et al. extended it to highorder methods using a continuous Chebyshev collocation method Rubio2013 () and later the Discontinuous Galerkin Spectral Element Method (DGSEM) Rubio2015 (). In that work, they studied the quasia priori truncation error estimation, which allows estimating the truncation error without having fully converged fine solutions, and introduced the concept of isolated truncation error (valid only for DG formulations), which only considers inner elemental contributions to the error and neglects the upwind contributions. More recently, Kompenhans et al. Kompenhans2016 () applied these estimators to perform panisotropic adaptation for the Euler and NavierStokes equations, and compared based to featured based adaption, showing better performance for the former Kompenhans2016a (). The adaptation strategy consisted in converging a high order representation (reference mesh) to a specified global residual and then performing a single error estimation followed by a corresponding padaptation process. Even though their methodology is very promising, we will show that it produces a large underestimation of the error for polynomial orders that are higher than the ones in the original reference mesh. This fact makes necessary to compute the initial solution in a very refined reference mesh to avoid inaccuracies.
In this paper, we extend the work on highorder estimators by Rubio et al. Rubio2013 (); Rubio2015 (); rubio2015truncation (), and formulate a new anisotropic truncation error estimator that exploits the tensor product basis expansion of the DGSEM. The new error estimator is shown to be suitable for performing anisotropic padaptation, and to have two main advantages over existing truncation error estimators; first, that it requires fewer operations to estimate the truncation error of all possible combinations of polynomial orders; and second, that it yields more accurate estimations of the truncation error for representations of a higher order than the reference mesh. This feature allows using reference meshes of a lower polynomial order, hence reducing the computational cost. We also analyze the properties of the traditional nonisolated truncation error and the isolated truncation error. To the authors’ knowledge, this is the first time that a highorder truncation error estimator based on the estimation method is formulated in an anisotropic/decoupled way, analyzed and tested.
The paper is organized as follows. In section 2, we present the mathematical background. First, the Discontinuous Galerkin Spectral Element Method is briefly summarized; then, we detail the existing techniques for approximating the truncation error of isotropic and anisotropic representations using the estimation method. In section 3, the proposed anisotropic estimator is introduced and analyzed. In section 4, we present a validation of the assumptions needed for formulating the new approach and study the properties of the proposed method by means of a manufactured solutions test case of the compressible NavierStokes equations. Finally, the conclusions are summarized in section 5.
2 Mathematical background
In section 2.1, we describe briefly the DGSEM approach. Section 2.2 contains the error definitions that will be used throughout the paper and provides an insight into the convergence properties of the different error measures. In section 2.3, we review the estimation method for DGSEM schemes, and then we explain in section 2.4 how it has been used in the literature for obtaining anisotropic error extrapolations.
2.1 The Discontinuous Galerkin Spectral Element Method (DGSEM)
We consider the approximation of systems of conservation laws,
(2) 
where is the vector of conserved variables, is the flux dyadic tensor which depends on , and is a source term. This system represents, among others, the compressible NavierStokes equations, as detailed in Appendix C. Multiplying equation 2 by a test function and integrating by parts over the domain yields the weak formulation:
(3) 
where is the normal unit vector on the boundary . Let the domain be approximated by a tessellation , a combination of finite elements of domain and boundary . Moreover, let , , and be approximated by piecewise polynomial functions (that are continuous in each element) defined in the space of functions
(4) 
where is the space of polynomials of degree at most defined in the domain of the element . Remark that the functions in may be discontinuous at element interfaces and that the polynomial order may be different from element to element. Equation 3 can then be rewritten for each element as:
(5) 
where the superindex “” refers to the functions as evaluated inside the element , i.e. ; whereas the superindex “” refers to the value of the functions on the external side of the interface . The numerical flux function, , allows to uniquely define the flux at the element interfaces and to weakly prescribe the boundary data as a function of the conserved variable on both sides of the boundary/interface ( and ) and the normal vector (). Multiple choices for the numerical flux functions can be found in the literature toro2013riemann (). In the present work, we use Roe Roe1981 () as the advective Riemann Solver and BassiRebay 1 Bassi1997 () as the diffusive Riemann solver. Remark that the numerical flux must be computed in a specific manner when the representation is nonconforming Kopriva2002 ().
Since , , and belong to the polynomial space , it is possible to express them inside every element as a linear combination of basis functions ,
(6) 
Therefore, equation 5 can be expressed in a discrete form as
(7) 
where is the local solution that contains the coefficients of the linear combination for the element ; is the global solution that contains the information of all elements; is known as the elemental mass matrix, and is a nonlinear spatial discrete operator on the element level:
(8)  
(9) 
Note that the operator is applied on the global solution, since it is the responsible for connecting the elements of the mesh (weakly). Assembling the contributions of all elements into the global system we obtain:
(10) 
In the DGSEM kopriva2009implementing (), the tesselation is performed with nonoverlapping hexahedral elements of order (independent in every direction) and the integrals are evaluated numerically by means of a Gaussian quadrature that is also of order . For complex geometries, it is most convenient to perform the numerical integration in a reference element and transform the results to the physical space by means of a highorder mapping:
(11) 
where the order of is at most (subparametric or at most isoparametric mapping). The differential operators can be expressed in the reference element in terms of the covariant () and contravariant () metric tensors:
(12) 
Under these mappings, the gradient and divergence operators become:
(13) 
where the Jacobian of the transformation can be expressed in terms of the covariant metric tensor:
(14) 
For details on how to compute the metric terms for 2D and 3D geometries, see Kopriva2006 (). Furthermore, in the DGSEM the polynomial basis functions ( in equation 2.1) are tensor product reconstructions of Lagrange interpolating polynomials on quadrature points in each of the Cartesian directions of the reference element:
(15) 
Therefore, are simply the nodal values of the solution, and is a diagonal matrix containing the quadrature weights and the mapping terms. In the present work, we make use of the LegendreGauss quadrature points kopriva2009implementing ().
2.2 Definition of Errors
In this section, we define some measures of the error that will be used throughout the paper.
Definition 1 (Interpolation error)
The difference between a function and its polynomial interpolant of order :
(16) 
where is the function that can be reconstructed using the polynomial expansion with coefficients (equation 2.1). For sufficiently smooth functions, in the asymptotic range the interpolation error in an element behaves as:
(17) 
where and are constants that depend on the local smoothness of the function canuto2012spectral (); hesthaven2007nodal () and is the local polynomial order in the element . In the DGSEM, the use of tensorproduct bases in dimensions allows decoupling the interpolation error in directional components, each of which depends solely on the polynomial order in the corresponding direction:
(18) 
As a consequence, in the panisotropic DGSEM, the interpolation error exhibits a tensorproducttype error bound in dimensions inside every element Rubio2015 (),
(19) 
Definition 2 (Discretization error)
The difference between the exact solution to the problem, , and the one obtained with a discretization of order , :
(20) 
The discretization error in an element is influenced by other elements because of the advection properties of the PDE. In fact, we will decouple the discretization error in locallygenerated and externallygenerated contributions for every element:
(21) 
In the pisotropic DGSEM, it can be assumed that the discretization error in each element behaves as Rubio2015 (); rubio2015truncation ():
(22) 
where is the number of elements, and are constants that depend on the smoothness of the solution in the element canuto2012spectral (); hesthaven2007nodal (), and and are constants that depend both on the smoothness of the solution and the advection properties of the PDE. The first term on the righthand side corresponds to the bound of the locallygenerated discretization error (), whereas the second term is the bound of the externallygenerated discretization error () which gathers the errors that are introduced through the Riemann solver. Note that is the minimum possible (lower bound of) . For an anisotropic representation in dimensions, the expression becomes Rubio2015 (); rubio2015truncation (); Georgoulis2003 ():
(23) 
Definition 3 (Quadrature error)
The quadrature error, also referred to as the numerical integration error, is the difference between the exact integral of a function and its approximation by a Gaussian quadrature:
(24) 
where the superindex on the integral indicates that it is approximated using a Gaussian quadrature of order ,
(25) 
and are the quadrature weights.
Definition 4 (Nonisolated truncation error)
We define the nonisolated truncation error of a discretization of order as the difference between the discrete partial differential operator of order and the exact partial differential operator applied to the exact solution:
(26) 
The exact partial differential operator can be derived from equation 2 as
(27) 
and the discrete partial differential operator is derived from equation 10 as
(28) 
where contains the sampled values of in all the nodes of the domain and is a sampling operator. is reconstructed from elementwise with equation 2.1. Since for steady state, , the nonisolated truncation error can be then computed inserting equation 28 into 26 as
(29) 
The dependence of the nonisolated truncation error on the discretization error is obtained by using definition 2 and expanding equation 29 as a Taylor series:
(30) 
Taking into account equations 30 and 23, and based on previous numerical results rubio2015truncation (); Rubio2015 (), Kompenhans et al. Kompenhans2016 () stated that the truncation error in an element is bounded by
(31) 
This expression was validated experimentally Kompenhans2016 (); Kompenhans2016a (). The first term in equation 31 is the bound of the locallygenerated truncation error, whereas the second term is the bound of the externallygenerated truncation error that enters through the Riemann solver and does not depend on the local polynomial orders. The second term is a consequence of the dependence of the discretization error on the solution in other elements.
Definition 5 (Isolated truncation error)
The isolated truncation error is defined as Rubio2015 ()
(32) 
where is the isolated discrete partial differential operator, which is derived in the same manner as , but is not substituted by during the process (equation 5). Therefore, the sampled form of the discrete isolated partial differential operator yields
(33) 
where the elemental contribution to the nonlinear discrete operator is
(34) 
This change eliminates the influence of the neighboring elements and boundaries in the truncation error of each element. The dependence of the isolated truncation error on the interpolation error of the fluxes inside an element () can be expressed as (see Appendix A and Rubio2015 ()):
(35) 
This shows that indeed depends only on the discrete representation of the numerical solution in the element . Rubio et al. Rubio2015 () pointed out that the isolated truncation error might be a better sensor for adaptation algorithms for hyperbolic PDEs than its nonisolated counterpart or the discretization error, since unlike the last two, the first one is not affected by neighbors’ errors. Notice that equation 35 resembles the DETE Roy2010 (). In this case, the isolated truncation error acts as a source term for the interpolation error. In consequence, decreasing the isolated truncation error reduces the interpolation error.
Finally, it is important to note that the spectral convergence of the isolated truncation error is similar to the nonisolated truncation error Rubio2015 (); rubio2015truncation () and can be expressed as
(36) 
Remark that, because of the reasons exposed above, in this case there is no externallygenerated truncation error. The hat notation will be dropped from now on since, unless explicitly stated, the formulations in this paper are valid for both the nonisolated and the isolated truncation errors.
2.3 Estimation method
Since in general the exact solution to the problem is not available, we are interested in using an estimation for equations 29 and 33. The estimation method makes use of an approximate solution on a reference mesh of order instead of the exact one. The most straightforward methodology is to converge this highorder solution to a low residual near machine roundoff, . This is known as the a posteriori approach. In practice, one can also use a nonconverged solution, . This is known as the quasia priori approach. In this paper, we use the following formulation, which is valid for the aposteriori method and the quasi apriori approach without correction:
(37) 
where is an interpolation operator from order to order . For compactness, the notation of this work omits the interpolation matrix such that . Equation 37 is valid for both the isolated (inserting the ) and the nonisolated truncation error. Note that the truncation error estimation can be easily performed for anisotropic polynomial representations of dimensions. For instance, in a 2D anisotropic case, equation 37 can be rewritten as:
(38) 
2.4 Low order extrapolation of the truncation error estimations
In this section, we review the method proposed by Kompenhans et al. Kompenhans2016 () to extrapolate the estimations of anisotropic representations. This method was successfully used to perform a padaptation strategy Kompenhans2016 (); Kompenhans2016a (). We will show that their strategy can be classified as a low order extrapolation. In order to do so, let us first introduce the concept of truncation error map.
Definition 6 (Truncation error map)
The (graphical) representation of the truncation error behavior inside an element with respect to the polynomial order as a dimensional plot of as a function of the polynomial order in every direction of the reference element , where is the number of dimensions.
Because of the spectral convergence of the truncation error in the asymptotic range, the onedimensional (or isotropic dimensional) map turns out to be a discrete scatter plot of points that describe a linear function with a negative slope (the convergence rate ), as shown in Figure 1(a).
Kompenhans et al. Kompenhans2016 () used the estimated truncation error map to adapt the polynomial orders of a given mesh using a specified maximum permitted error threshold, . The method for estimating the map consists of four steps:

Use the inner map to look for a combination of polynomial orders that fulfills the specified error threshold. If a combination fulfills , adapt the polynomial order and exit the adaptation process. Otherwise, additional considerations are required.

Compute and perform a linear regression analysis in the direction in order to describe the behavior of as a function of (). The result of the linear regression is marked with a dashed line in Figure 1.

Use the linear regression to estimate the truncation error for , and select the value of and independently from these extrapolations. The extrapolated values of the truncation error are marked with red squares in Figure 1.
This procedure is performed for every element in all the Cartesian directions to adapt the mesh. For further details, refer to the original paper by Kompenhans et al. Kompenhans2016 () and to our example in section 4.2.
2.4.1 Analysis of the method
Two remarks can be made about the described fourstep procedure:
Remark 1
Steps 3 and 4 assume that the spectral convergence observed in 1D extends to higher dimensions along iso lines of the truncation error map.
Remark 2
For the nonisolated truncation error, the behavior shown in Figure 1 can only be expected for the locallygenerated component (see equation 31). This means that the extrapolation procedure may predict unexpected behaviors if the truncation error in neighboring elements is high, i.e., if the estimation procedure is not performed elementwise while keeping the polynomial order in other elements sufficiently^{1}^{1}1Sufficiently high does not necessarily mean that the polynomial order of the other elements must be kept in , but that it must be high enough so that the externallygenerated contributions to the truncation error are less than the internallygenerated ones. high.
As stated in the remark 1, the extrapolation procedure assumes that the truncation error decreases exponentially along iso lines of the truncation error map. That is the same as saying that the truncation error map is a plane for , and in general that it is a hyperplane of dimension . In Figure 2 we present an illustration that resembles the hyperplane behavior in two dimensions for perfect spectral convergence. The described methodology consists in constructing iso lines on the hyperplane, which should contain the values of the truncation error for (red line with triangular markers and black line with circular markers in Figure 2). In that scenario, selecting independently can be regarded as a conservative criterion, since in the hyperplane we have:

In 2D:

In 3D:
for . See Figure 2.
Hereinafter, the method by Kompenhans et al. shall be referred to as the low order extrapolation method, since it supposes that the truncation error map () has a linear behavior. In light of the analysis in section 3, we will be able to formulate a high order extrapolation method that provides extrapolated estimations with increased accuracy.
3 New anisotropic truncation error estimation
In this section, we present the new anisotropic truncation error estimator, discuss some of its properties and compare them with the error estimators that have been used in the literature for performing anisotropic padaptation. The formulation of this new estimator involved a mathematical proof based on specific assumptions that is detailed in section 3.1. In section 3.2, we analyze the convergence behavior of the anisotropic estimator. In section 3.3, we describe how the new estimator can be used for approximating the truncation error of higherorder representations, and prove that it is superior to existing estimators at obtaining these approximations.
3.1 Anisotropic estimation
The anisotropic estimator is a generalization of the ideas reviewed in section 2.3 and is based on four assumptions that are explained first. For the sake of readability and without loss of generality, all the mathematical statements in this section are for 2D formulations, where are the polynomial orders in the 2 directions of the reference element. However, all the statements and proofs can be directly generalized to dimensions.
Assumptions
Following assumptions are a consequence of the tensor product basis functions of the DGSEM and hold for sufficiently smooth solutions in the asymptotic range. The assumptions are:

The truncation error has an anisotropic behavior and, therefore, can be decoupled in its directional components:
(40) Here, it is important to note that is the projection of the global truncation error, , into the local direction, .

The locallygenerated truncation error in each direction depends only on the polynomial order in that direction:
(41)
Assumptions (a) and (b) follow from the work of Rubio et al. Rubio2015 (); rubio2015truncation (). Furthermore, assumption (b) relates to the anisotropic spectral convergence behavior of the truncation error (equations 31 and 36).
Theorem 1
The truncation error of a DGSEM discretization of order can be approximated from a semiconverged solution of order , such that , as the sum of the directional estimations obtained by coarsening in the different space dimensions:
(42) 
Proof
This proof is specific for the isolated truncation error. We refer to Appendix B for a brief proof that Theorem 1 also holds for the nonisolated truncation error under additional assumptions.
Let us note that assumptions (a) and (b) are consistent with the dependence of the isolated truncation error on the interpolation error (equation 35) and the anisotropic behavior of the latter (equation 19).
We start by obtaining the analytical expression for the isolated estimation. To that end, we use the same procedure as in Appendix A. The estimate of the isolated truncation error in the DGSEM can be expressed for any basis function in an element as
(43) 
In this case, instead of the exact solution to the problem, , we use a solution on a higher order mesh, . Therefore, using the definition of interpolation error (equation 16) and discretization error (equation 20), the flux is
(44) 
and the source term is
(45) 
(46) 
Remark that although the isolated truncation error of an element does not depend on external sources, its approximation by estimation is affected by the discretization error of the reference mesh, . This translates into a weak influence of external (upwind) errors transmitted through the Riemann solver. In the twodimensional case and coarsening in only one direction (here the direction ), equation 46 becomes
(47) 
Now, we rewrite equation 47 decoupling the interpolation error in directional components and taking into account that (equation 18):
(48) 
Notice that all terms on the righthand side, except for the first one and the quadrature error, are of the order of errors on the higherorder mesh. Therefore, and taking into account that we are coarsening in the direction , for sufficiently smooth solutions we can expect the first term on the righthand side to be the leading term. Simplifying, the directional estimation provides
(49) 
On the other hand, inserting equation 18 into 35, the isolated truncation error of a representation of order yields
(50) 
Notice, again, that the directional components of the interpolation error only depend on the polynomial order in the corresponding direction (equation 18). Therefore, neglecting additional quadrature errors, we recover equation 42 for the isolated truncation error by combining equations 49 and 50:
(51) 
From the previous analysis we can conclude that, when the estimation method is performed coarsening only in the direction , the result is an approximation to the directional component of the truncation error, . We can also arrive at this conclusion intuitively if we realize that cannot be better than at describing the solution in the direction .
Theorem 1 can be easily generalized to three dimensions to obtain
3.2 Convergence behavior of the anisotropic truncation error
In this section we analyze the convergence properties of the truncation error map using Theorem 1. Let us first consider the directional components of the truncation error.
Theorem 2
The directional components of the locallygenerated truncation error exhibit spectral convergence with respect to the polynomial order in the corresponding direction:
(52) 
Proof
According to assumption (equation 41), each directional component of the locallygenerated truncation error, , depends solely on the polynomial order in the corresponding direction, . If we insert equation 41 into equation 36 and analyze the dependencies term by term, we recover 52 for the isolated truncation error. In the same way, if we insert equation 41 into 31, we recover 52 for the nonisolated truncation error.
Now, we are able to analyze the convergence behavior along lines of the truncation error map, where a line in dimensions is defined as:
(53) 
with .
Theorem 3
The total truncation error does not necessarily decrease exponentially along lines of the truncation error map.
Proof
The corresponding positive statement can be easily proven wrong with a counterexample. Let us consider a 2D anisotropic representation. From theorem 2, we know that the truncation error of each directional component in an element, , decreases exponentially when increasing ; and that the decreasing rate is , a constant that depends on the smoothness of the solution in the direction . Let us suppose that for a certain element in a mesh, the directional components of the error have the same value for a specific combination of polynomial orders in a certain norm:
(54) 
Let us analyze the convergence rate along an iso line of the truncation error map with constant , i.e. the convergence rate of as a function of . Note that, according to assumption (b), along the line we have
(55) 
Furthermore, assumption (a) states that the total truncation error along the line of the map is:
(56) 
(57) 
Remember that the truncation error map is defined as the dependence of on (definition 6). Therefore, taking logarithms in both sides and rearranging, equation 57 can be rewritten in two equivalent forms:
(58)  
(59) 
For , the second term on the righthand side of equation 58 vanishes, which indicates that the convergence rate of the truncation error along an iso line of constant tends to . On the other hand, for , the second term on the righthand side of equation 59 vanishes, which implies that the convergence rate along an iso line of constant tends to zero, since the truncation error is bounded by , a constant that does not depend on . In other words, the iso line on the hyperplane for is not a straight line.
3.3 High order extrapolation of the truncation error estimations
In this section, we present a procedure for extrapolating the truncation error estimations (inner map) that can be obtained by applying Theorem 1. Since Theorem 3 rules out the possibility of extrapolating along iso lines of the truncation error map, we take advantage of the anisotropic behavior of the truncation error (equation 40) and the spectral convergence of its directional components (Theorem 2). The proposed methodology, which is valid for both the isolated and nonisolated truncation error, can be summarized in three steps:

Perform anisotropic coarsening to obtain (Theorem 1) and construct the inner truncation error map directly. In dimensions, this requires only evaluations of the discrete partial differential operator , where
(60) 
Compute and perform a linear regression analysis in the direction in order to describe the behavior of as a function of . This is supported on the proved spectral behavior of the directional components of the truncation error (Theorem 2).

Construct the outer truncation error map using equation 42 and the extrapolated values of :
(61)
An example is provided in section 4.
3.4 Theoretical comparison of the new anisotropic estimation with previous approaches
Figure 3(a) illustrates the theoretically predicted behavior of the truncation error map that is obtained with the new anisotropic estimation method for a toy problem (only illustrative). As noted in the proof of Theorem 3, remark that along an iso line, the truncation error firstly decays exponentially for low and then tends asymptotically to a constant value for high . As can be seen, contrary to the loworder extrapolation, the new anisotropic estimation does not assume that the truncation error map has a linear behavior. That is the reason why it is called highorder extrapolation. Figure 3(b) shows a comparison of the hyperplane behavior with the one obtained using the new anisotropic estimation method. As can be seen, the hyperplane tends to underpredict the truncation error for some combinations of polynomial orders as compared to the new truncation error estimator. A comparison of the output of both estimation methods with the exact truncation error in a test case is provided in section 4.2.
It is noteworthy that even though the method of Kompenhans et al. Kompenhans2016 () supposes hyperplane behavior, their strategy of selecting the polynomial order in every direction independently minimizes the error involved in the low order extrapolation (compare the values of in Figures 3(a) and 2). Furthermore, remark that for and the new anisotropic estimation tends to have a hyperplane behavior. In fact, Kompenhans et al. (Kompenhans2016, , section 5.1) state that only the values of the truncation error where or should be used for the least square fitting. The highorder extrapolation can be seen as a form of bypassing this requirement.
Table 1 provides a final summary comparing the new estimation method and the previous methodology by Kompenhans et al. Kompenhans2016 ().
Feature  Kompenhans et al. Kompenhans2016 ()  Proposed estimation 

Number of evaluations of for inner map  
Accuracy of inner map  Very good  Good 
Accuracy of outer map  Poor  Good 
4 Validation of the anisotropic estimation method
The compressible NavierStokes equations can be written in conservative form (see Appendix C) and discretized using the DGSEM, as explained in section 2.1. In order to test the accuracy of the proposed estimation method, a 2D manufactured solutions test case is analyzed. The exact solution selected for the problem is
(62) 
which is simulated in the unit square, as depicted in Figure 4. Inserting equation 4 into 81, the source term for the 2D compressible NavierStokes equations yields
(63) 
The main interest is to validate the proposed error estimator and compare its outcome with previous works. Since the method of Kompenhans et al. Kompenhans2016 () explained in section 2.4 was formulated and used with the nonisolated truncation error, the results that are shown in sections 4.1 and 4.2 were obtained for the nonisolated truncation error estimator (). However, similar results can be obtained for the isolated truncation error, since its maps exhibit the same behavior as the ones presented here. In addition, in section 4.3 we compare the truncation error estimator with the isolated truncation error estimator when used for driving a padaptation procedure.
4.1 Truncation error maps and number of degrees of freedom
A fully timeconverged solution () of order 5 () is used to estimate the truncation error using the method of section 3.1. The results for element A are depicted in Figure 5(a). Figure 5(b) shows the exact truncation error. It can be seen that the proposed method predicts a truncation error map that is very similar to the exact one, even for extrapolated values. Hence, in agreement with the obtained results, the assumptions of section 3.1 are reasonable. These maps can be used for selecting an appropriate combination of polynomial orders such that a maximum truncation error is achieved employing a minimum number of degrees of freedom.
Figure 6 shows the map of the number degrees of freedom (DOFs) for every (,)combination. The polynomial orders that achieve a truncation error are marked with black squares. Let us remark that, although these results are not exactly the same for the estimated and exact truncation error maps, they are very similar. Therefore, we conclude that the proposed estimation method may be used for adaptation purposes. Notice that there are many alternatives that produce a truncation error in the desired range, but there is only one that minimizes the number of degrees of freedom and, therefore, the computational cost.
4.2 Comparison with previous methodologies
Figure 7 shows the 3D representation of the exact truncation error map (a), the one obtained with the highorder extrapolation (b), and the one obtained with the loworder extrapolation (c) here, we illustrate the complete hyperplane. The maps were generated with the same fully timeconverged solution of order . As can be seen, the truncation error map generated with the highorder extrapolation bears close resemblance to the exact one, whereas the hyperplane underpredicts the truncation error in some regions, as anticipated in section 3.4.
If we generate the truncation error map using the method of Kompenhans et al. Kompenhans2016 () (section 2.4), we obtain Figure 10. Remark that although the method of Kompenhans et al. produces accurate results for , it fails to predict the behavior of the truncation error for . In fact, using this method the full truncation error map is not being generated, but only the extrapolations for the iso lines and .