(traditionally in hexahedra) and stores the data at the Gauss or Gauss-Lobatto nodes of the quadrature
rule. The use of a quadrature rule with the same number of nodes as the approximate solution equips the
DGSEM with a diagonal mass matrix and very cheap-to-compute operators. In fact, the computational
cost of the DGSEM has been estimated to be a factor of four smaller than other conventional DG methods
[9]. In addition, since the DGSEM uses tensor-product bases, it can handle p-anisotropic discretizations
efficiently [10,11], i.e. discretizations that have different polynomial degrees in each coordinate direction.
For all those properties, the DGSEM has been used in a wide range of applications, including the simulation
of incompressible Navier-Stokes, compressible Navier-Stokes, Cahn Hilliard equation, or multi-phase flows
(see [12] and references therein).
In their famous review paper, Wang et al. [3] point out that one of the challenges the high-order commu-
nity must address to impact the design process and replace traditional low-order codes is the development
of efficient mesh adaptation strategies. The idea behind these strategies is to reduce the number of degrees
of freedom (DOFs) while maintaining high accuracy, which translates into shorter computational times and
reduced storage requirements. Local adaptation can be performed by subdividing or merging elements (h-
adaptation), by enriching or reducing the polynomial degree in certain elements (p-adaptation), by relocating
the position of the nodes in a mesh (r-adaptation). For all these strategies it is of paramount importance
to identify the flow regions that require refinement or coarsening with a local error estimation.
Adaptation strategies have been classified according to the type of error measure that is employed as
feature-based adaptation, adjoint-based adaptation, and local error-based adaptation. A comparison of these
three approaches was performed by Fraysse et al. [13] for finite volume approximations and by Kompenhans
et al. [14] and Naddei et al. [15] for high-order DG methods. The feature-based adaptation is the classical
approach and uses easy-to-compute error measures that depend on the flow features. They rely on the
assumption that high errors are expected where the flow is more difficult to resolve. Hence, refinement is
predicted where high velocity, density or pressure gradients are identified [16,17]. For DG discretizations,
an easy-to-compute feature-based adaptation criterion is the assessment of jumps across element interfaces
[18–20]. The main disadvantage of these methods is that there is no direct relation between the adaptation
criterion and the numerical errors and thus the accuracy is not easily predictable. Additionally, the only
way to solve steady-state problems is to adapt iteratively.
A second and more sophisticated approach is known as adjoint-based adaptation. In this approach, a
functional target is defined (e.g. drag or lift in external flow aerodynamics) and the adjoint problem is solved
to obtain a spatial distribution of the functional error, which is then used to adapt the mesh. This technique
was originally developed for structural analysis using FEM by Babuˇska and Miller [21,22], and has been
used recently for adaptation strategies in DG methods [23–25]. The main drawback of this approach is the
high computational cost to solve the adjoint problem and the storage requirements needed to save the error
estimators, especially in unsteady flows. Moreover, only the error of the functional analyzed is guaranteed
to be reduced, whereas the error of other functionals may deteriorate.
A computationally more efficient alternative is the local error-based adaptation, which is based on the
assessment of any measurable (not feature-based) local error in all the cells of the domain [24]. The local
error-based adaptation methods are interesting since, in contrast to feature-based methods, they provide a
way to predict and control the overall accuracy, and are computationally cheaper than adjoint-based schemes
[10,14]. For those reasons, local error-based adaptation strategies are retained in this work. A large amount
of effort has been invested in the development of reliable local error-based adaptation methods. Estimations
of the local discretization error have been used by Mavriplis [26,27] to develop hp-adaptation techniques for
the spectral element method. Residual-based p-adaptation is also a local error-based adaptation method,
which uses the residual to measure how accurate is the local approximation. This method was originally
developed for Finite Elements (FE) and has been successfully used with DG methods [15,23]. In the case
of modal (hierarchical) DG methods, a possibility is to employ low cost error estimates that take advantage
of the modal approximation to drive p-adaptation procedures, such as the Variational Multiscale (VMS)
indicator by Kuru and De la Llave Plata [28], or the spectral decay indicator by Persson and Peraire [17].
In this work, we favor truncation error estimators, another local error-based alternative to drive a mesh
adaptation method. The truncation error is related to the discretization error through the Discretization
Error Transport Equation [29], where it acts as a local source term. This relation makes it useful as an
2