Javier García Bautista
email: javier.garcia_bautista@onera.fr.

Marta De La Llave Plata
email: marta.de_la_llave_plata@onera.fr.

Vincent Couaillier
email: vincent.couaillier@onera.fr.

Michel Visonneau
email: michel.visonneau@ec-nantes.fr.

Kai Schneider
email: kschneid@univ-amu.fr.

Adaptive Discontinuous Galerkin Approach Based on Local Multiwavelets Analysis

The main goal of the proposed research is to investigate and develop error estimators and mesh adaptation strategies within a discontinuous Galerkin (DG) formulation in order to reduce the computational cost, for a prescribed level of accuracy. We are interested in how multiwavelets (MWs) and their properties may shed new light on the adaptation process. This is motivated by the fact that MWs break any input apart into a hierarchy of low resolution data and subsequently finer details. Our error estimator makes use of MWs' properties while being local to the element, maintaining the parallel efficiency of the solver. Early tests on the one-dimensional viscous Burgers' equation show promising results. This work will be focused on a backward-facing step configuration to assess the performance of the method in higher dimensions.

I. Introduction

I the framework of the development of the CFD solver Aghora [1-3], ONERA is working on the development of hp-adaptive high-order discontinuous Galerkin (DG) methods.

Finite volume methods (FVM) have been extensively studied and developed by the CFD community. The FVM uses piecewise constant basis functions defined within the element. On the other hand, finite element methods (FEM) employ basis functions of different orders with continuous global support. DG approaches appear as a hybrid alternative to the aforementioned methods: the basis functions are local to each element while providing high-order accuracy.

Moreover, the DG method offers the possibility to efficiently adapt the spatial resolution by either modifying the local mesh size (h-adaptation), the local polynomial degree (p-adaptation), or both simultaneously (hp-adaptation). In h-adaptation strategies, flagged elements are divided into smaller elements (h-refinement) or agglomerated into a larger element (h-coarsening). On the other hand, p-adaptation approaches increase (p-enrichment) or lower (p-coarsening) the degree of the polynomial approximation in marked elements. An hp-decision algorithm selects the most appropriate adaptation strategy (h-and/or p-) to be adopted for each element. Smooth regions in the solution promote p-adaptation, whereas regions featuring locally steep gradients (e.g. shocks and boundary layers) are better captured using h-adaptation. Overall, the process should maintain the same accuracy than the uniformly refined grid (for the same minimum mesh size), yet for a lower number of degrees of freedom (DOFs).

With this in mind, the research in [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Discontinuous Galerkin Methods for the Euler and Navier-Stokes Equations[END_REF][START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF] shows the benefits of p-adaptation, and how its use can considerably reduce the number of DOFs for a prescribed level of accuracy, compared to the use of a uniform approximation order.

In this work the focus is on h-adaptation strategies, and in particular on the formulation and implementation of multiwavelet-based (MW) error estimators (EE). Multiwavelets belong to the broader field of multiresolution analysis (MRA) [START_REF] Harten | Multiresolution Representation of Data: A General Framework[END_REF]. This theory refers to the possibility of representing some given data into scale-wise contributions. The aggregation of the coarse scale and the nested sequence of details yields the original data themselves. This scale separation represents the main idea behind wavelets (scalar MRA theory) and multiwavelets (vectorial MRA theory).

In the field of CFD, the natural step forward is for the MRA to represent the numerical solution of conservation laws. On a grid it is easy to visualise how more and more details of the conserved quantities can be extracted by refining the mesh until the flow field is fully resolved. The MRA should be adapted to the numerical method in use.

A fully adaptive FVM approach based on MRA has been developed in [START_REF] Müller | Adaptive Multiscale Schemes for Conservation Laws[END_REF][START_REF] Müller | Multiresolution Schemes for Conservation laws[END_REF] by explicitly using biorthogonal wavelets to represent the numerical solution. Another line of research, consists in employing the recursive use of projection and prediction operators to define the coarse and fine levels of the cell-averaged values [START_REF] Roussel | A Conservative Fully Adaptive Multiresolution Algorithm for Parabolic PDEs[END_REF][START_REF] Domingues | Space-Time Adaptive Multiresolution Methods for Hyperbolic Conservation Laws: Applications to Compressible Euler Equations[END_REF][START_REF] Deiterding | Adaptive Multiresolution or Adaptive Mesh Refinement? A Case Study for 2D Euler Equations[END_REF][START_REF] Deiterding | Comparison of Adaptive Multiresolution and Adaptive Mesh Refinement Applied to Simulations of the Compressible Euler Equations[END_REF]. These are all examples of the use of scalar MRA theory in combination with FVM to perform mesh adaptation.

On the other hand, when dealing with DG the vectorial MRA theory comes as the natural choice. Multiwavelets can be easily combined with DG due to its flexibility in matching the order of the approximation, while keeping compact support [START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF].

Following these ideas, MWs and DG have been successfully combined to perform mesh adaptation in [START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF][START_REF] Hovhannisyan | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws[END_REF][START_REF] Gerhard | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws: multidimensional case[END_REF][START_REF] Gerhard | Multiwavelet-Based Grid Adaptation with Discontinuous Galerkin Schemes for Shallow Water Equations[END_REF]. In the approach proposed by the authors, the entire physical domain is subjected to the MW decomposition. That is, each level in the MW space corresponds to a refinement level in the physical space. On a different line of research, the work of [START_REF] Vuik | Multiwavelet Troubled-Cell Indicator for Discontinuity Detection of Discontinuous Galerkin Schemes[END_REF][START_REF] Vuik | Automated Parameters for Troubled-Cell Indicators Using Outlier Detection[END_REF] combines DG and MWs for the detection of singularities. The most extensively used MWs are the so called Alpert's multiwavelets [START_REF] Alpert | A Class of Bases in 2 for the Sparse Representation of Integral Operators[END_REF][START_REF] Alpert | Adaptive Solution of Partial Differential Equations in Multiwavelet Bases[END_REF], in which the MW bases are associated to scaled Legendre polynomials. We will also employ this form of MW in our method, which will be presented in the next sections.

As discussed in the previous paragraph, in the DG h-adaptation literature, references [START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF][START_REF] Hovhannisyan | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws[END_REF][START_REF] Gerhard | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws: multidimensional case[END_REF][START_REF] Gerhard | Multiwavelet-Based Grid Adaptation with Discontinuous Galerkin Schemes for Shallow Water Equations[END_REF] constitute a representative example of h-adaptation based on global-multiwavelets. Namely, a cascade of scales is built from the DG solution on the finest mesh (reference mesh). This new representation sheds light on the local structure of the solution. Indeed, the detail coefficients become small with incrementing resolution level when the analysed solution is smooth. Consequently, thresholding may be used to perform local grid adaptation in those elements in which the MW coefficients still maintain a significant value (refining), and detect the regions where the regularity of the solution allows us to decrease the local resolution while maintaining the desired level of accuracy (coarsening). Pursuit of this philosophy produces the final adapted mesh (compression of DG solution).

This method, while soundly based on wavelet theory and accurate, bears important constraints. Firstly, a solution on a finer mesh must be foreknown. Secondly, Cartesian grids in which a dyadic hierarchy can be built must be employed. Lastly, parallelisation may prove challenging due to the global character of the MRA technique.

To overcome these shortcomings, in the present work MW analysis is performed locally within each element. By being local to the element, adaptation can be applied by starting from a coarse mesh. More general grids may also be used (not limited by the strict translation and dilation properties of MWs) and the parallel efficiency of the original DG method is conserved.

However, the original DG element-wise solution must be enriched with new data so that the local multiwavelet expansion may extract significant information. This is done by employing a reconstruction process involving the current element and its neighbours. The MRA decomposition is then locally performed on the post-enriched DG solution. The resulting multiwavelet coefficients constitute the backbone of an error estimator which will drive the adaptation process.

Therefore, the main goal of this work is to present a new error estimator grounded on an local multiwavelet expansion of a DG solution which is exposed to enrichment by a post-processing reconstruction. The error estimator is tested for the one-dimensional viscous Burgers' equation and a two-dimensional laminar flow over the backward-facing step.

In section II the DG method is explained for the 1D Burgers' equation and the 2D compressible Navier-Stokes equations. In the next section the theoretical background on multiwavelets in one and two dimensions is presented. Moreover, the connection between DG and MWs is established. Section IV shows the structure of the proposed MW-based error estimator. The enrichment of the original DG solution and the construction of the local MRA hierarchy are also described in this section, together with the description on how the error estimator is evaluated and of the element marking strategy. Finally, in section V numerical results are discussed for the 1D Burgers' equation and the 2D laminar backward-facing step flow.

II. Modal Discontinuous Galerkin

In our work the one-dimensional unsteady viscous Burger's equation and the two-dimensional compressible Navier-Stokes equations will be studied. Let Ω ⊂ R be a bounded domain, where is the spatial dimension. Given appropriate boundary conditions on Ω, they can be written under the general expression

u + ∇ • (F (u) -F (u, ∇u)) = 0 , ∀x ∈ Ω, > 0 , (1)
u(x, 0) = u 0 (x), ∀x ∈ Ω , (2)
where u is the vector of conservative variables. The vectors F , and F are the convective and viscous fluxes, respectively. The DG method is based on the discrete weak formulation of the problem given in Eqs. [START_REF] Chapelier | Evaluation of a High-Order Discontinuous Galerkin Method for the DNS of Turbulent Flows[END_REF][START_REF] Chapelier | Inviscid and Viscous Simulations of the Taylor-Green Vortex Flow Using a Modal Discontinuous Galerkin Approach[END_REF]. We now proceed to partition the domain Ω into a shape-regular grid, Ω ℎ , formed by non-overlapping and non-empty elements of characteristic size ℎ. Interior and boundary faces in Ω ℎ are defined by E and E , respectively, such that E ℎ = E ∪ E . We search for approximate solutions in the function space of piecewise polynomials

V ℎ = { ∈ 2 (Ω ℎ ) : | ∈ P ( ), ∀ ∈ Ω ℎ } with degree at most . Let ( 1 , . . . ,
) ∈ P ( ) be a hierarchical and orthonormal modal basis of V ℎ , with dimension . The element-wise solution is then expressed as

u ℎ (x, ) = ℓ=1 U ℓ ( ) ℓ (x), ∀x ∈ , ∈ Ω ℎ , ∀ > 0 , (3)
where the coefficients (U ℓ ) 1≤ℓ ≤ are the degrees of freedom (DOFs) representing the solution in element . In this work we will only consider Cartesian meshes and thus the basis, ( ℓ ) 1≤ℓ ≤ , corresponds to the Legendre polynomials.

A mapping is then established between physical and reference space, Ξ = [-1, 1] . The variational form of Eq. ( 1) then follows: find

u ℎ in V ℎ such that ∀ ℎ ∈ V ℎ , for each element , we have ∫ u ℎ ℎ dΩ + L (u ℎ , ℎ ) + L (u ℎ , ℎ ) = 0 . (4)
The convective and viscous terms, L and L respectively, will be described in the following sections. To that end, additional notation is required: for a given interface in E , we define the average operator u = (u + + u -)/2 and the jump operator [[u]] = (u +u -) ⊗ n, where u + and u -are the traces of u at the interface between elements + and -, and n represents the orthonormal outward direction to an element + .

A. Viscous Burgers' Equation

In this case Ω ⊂ R and the vector of conservative variables simply becomes u = . The convective and viscous fluxes take the form F ( ) = 1 2 2 and F ( , ) = , respectively. Finally, the variational projection of the convective and viscous terms in Eq. (4) yields

L ( ℎ , ℎ ) = [ℎ ℎ ] - ∫ 1 2 2 ℎ ℎ d , (5)
L ( ℎ , ℎ ) = -[Θ ℎ ] - ∫ ℎ ℎ d -(ℎ -ℎ ) ℎ , (6)
in which the structure of presented in [START_REF] Alhawwary | On the Accuracy and Stability of Various DG Formulations for Diffusion[END_REF] has been followed. The approximation of the convective flux on an element face, ℎ , is fully defined by the local Lax-Friedrichs flux (LLF) [START_REF] Cockburn | An Introduction to the Discontinuous Galerkin Method for Convection-Dominated Problems[END_REF], i.e.

ℎ = F ( ) - 1 2 LLF [[ ]], LLF = max min( -, + ) ≤ ≤max( -, + ) dF ( ) d . ( 7
)
The numerical viscous flux, given by ℎ and Θ , is approximated by the symmetric interior penalty method [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF]. It reads

ℎ = ℎ , (8)
Θ = ℎ -SIP [[ ℎ ]] . (9)
The penalty parameter, SIP , depends on the size of the element and the polynomial degree .

For the time integration an explicit time marching scheme is employed, the strong stability preserving (SSP) 3rd-order 4-stages Runge-Kutta method [START_REF] Carpenter | Fourth-Order 2N-storage Runge-Kutta Schemes[END_REF]. Consequently, Eq. ( 4) can be expressed as

ℓℓ ℓ = L ℓ ℓ -1 , ℓ , ℓ +1 , ℓℓ = ∫ ℓ ℓ d , (10)
where L ℓ ℓ -1 , ℓ , ℓ +1 encompasses the convective and viscous terms of Eqs. [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF][START_REF] Harten | Multiresolution Representation of Data: A General Framework[END_REF], and the DOFs are the unknowns of the system of ordinary differential equations.

B. Compressible Navier-Stokes Equations

We define Ω ⊂ R 2 with u = ( , v, ) the vector of conservative variables. The velocity vector is given by v = ( 1 , 2 ) , and = ( -1) + v•v 2 represents the specific total energy. The static pressure, , is defined by the ideal gas law, = , and = > 1 is the ratio of specific heats. Convective and viscous fluxes become

F = , v ⊗ v + Ī, ( + )v , (11)
F = 0, ¯ , ¯ • v -q , (12)
where ¯ = ∇v + (∇v) -2 3 (∇ • v) Ī represents the shear-stress tensor, in which is the dynamic viscosity as defined by the Sutherland's law. Lastly, the quantity q = -∇ is the heat-flux vector, with being the thermal diffusivity.

The discrete variational form of the convective terms in Eq. ( 4) reads

L (u ℎ , ℎ ) = ∫ F (u ℎ ) • ∇ ℎ dΩ + ∫ E [[ ℎ ]]h (u + ℎ , u - ℎ , n) dE + ∫ E + ℎ F u • n dE , (13)
where the boundary values u = u u + ℎ , u ext , n , with u ext a reference external state, are computed so that the boundary conditions are satisfied on E . The convective flux on an element face is approximated by the numerical flux h and it must satisfy the conditions of consistency and conservativity [START_REF] Cockburn | An Introduction to the Discontinuous Galerkin Method for Convection-Dominated Problems[END_REF]. In this work, the LLF flux is employed:

h u + ℎ , u - ℎ , n = F (u ℎ ) • n + 1 2 LLF u + ℎ -u - ℎ , LLF = max (J ( )) : = ± ℎ , (14)
where J ( ) = ∇ (F ( ) • ) indicates the Jacobian matrix of the convective fluxes in the direction and its spectral radius.

Finally, for the discrete variational form of the viscous terms, the BR2 method presented in [START_REF] Bassi | A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations[END_REF] is employed:

L (u ℎ , ℎ ) = ∫ F (u ℎ , ∇u ℎ + L ℎ ) • ∇ ℎ dΩ - ∫ E [[ ℎ ]] F u ℎ , ∇u ℎ + BR2 l ℎ • n dE + ∫ E + ℎ F u , ∇u + BR2 l ℎ • n dE , (15)
where the boundary values u , ∇u are consistent with the boundary conditions imposed on E and BR2 is a user-defined parameter for the stabilization of the method. The local lifting operator, l ℎ , extends its support on the current element and its neighbours, that is supp(l ℎ ) = { + ∪ -}. It is defined for each internal faces = E , = + ∩ -, as follows:

∫ + ∪ - ℎ l ℎ dΩ = - ∫ [[ ℎ ]] [[u ℎ ]] dE . ( 16
)
An analogous equation consistent with the boundary conditions can be obtained for the boundary faces E . The global lifting operator L ℎ is defined for the element as the sum of the local lifting operators, L ℎ = ∈ l ℎ . Finally, once every term has been defined, Eq. ( 4) results in a non-linear system of ordinary differential equations in which the polynomial coefficients, U ( ), are the unknowns to be determined. In this paper the work is focused on steady problems. Unlike the one-dimensional setting, in which an explicit scheme has been used, to evolve the solution in time we use the implicit Euler method instead. The subsequent non-linear system of equations is solved by the GMRES method with ILU(0) preconditioning.

III. Multiwavelets

In this section the basics of multiwavelet theory will be presented. We will proceed in a similar manner to previous section: first the concepts will be described in a one-dimensional framework and then we will move to higher spatial dimensions.

MRA and multiwavelets

Multiresolution analysis allow us to decompose a given input signal into a hierarchy of approximations of that signal at different levels of resolution [START_REF] Harten | Multiresolution Representation of Data: A General Framework[END_REF]. The level number will be denoted by . How much detail is captured by a particular resolution level depends on how many subdivisions, , this level owns. If we suppose that the relation between level and number of subdivisions is dyadic and we work in 2 ( [-1, 1]) then the support of the subdivisions is determined by

= [-1 + 2 -+1 , -1 + 2 -+1 ( + 1)], = 0, ..., 2 -1. (17)
Now we define two subspaces connected to each , denoted by , . The orthonormal bases for each subspace are called scaling functions, V = span{ ℓ }, and multiwavelets, W = span{ ℓ }, with ℓ = 0, ..., . This is the premise behind Alpert's multiwavelets [START_REF] Alpert | A Class of Bases in 2 for the Sparse Representation of Integral Operators[END_REF][START_REF] Alpert | Adaptive Solution of Partial Differential Equations in Multiwavelet Bases[END_REF]. The following development is based on this family of multiwavelets.

We start from the coarse level = 0 and build up from there. In this case the subspace of scaling functions results in V 0 ( [-1, 1]) and its basis [START_REF] Alpert | Adaptive Solution of Partial Differential Equations in Multiwavelet Bases[END_REF]:

0 ℓ0 ( ) = ℓ ( ) = (2ℓ+1) 2 ℓ ( ), ∈ [-1, 1] 0, otherwise (18)
where ℓ ( ) indicates the Legendre polynomials. Scaling functions in Eq. ( 18) up to = 2 are plotted in Fig. [START_REF] Chapelier | Evaluation of a High-Order Discontinuous Galerkin Method for the DNS of Turbulent Flows[END_REF]. On the other hand, multiwavelets undergo a more complex building process. The algorithm starts with a piecewise monomial of degree ℓ = 0, ..., in [-1, 1]. A Grand-Schmidt orthonormalisation is followed by an operation to increase the number of vanishing moments of the resulting function. This is enforced by ensuring orthogonality with respect to a higher degree monomial. The complete algorithm can be found in [START_REF] Alpert | A Class of Bases in 2 for the Sparse Representation of Integral Operators[END_REF]. The multiwavelets that spans the subspace W 0 ( [-1, 1]) are thus formed by these orthonormal functions ℓ, ( ) as follows:

0 ℓ0 ( ) = ℓ ( ) =          (-1) ℓ+ +1 ℓ, (-), ∈ [-1, 0] ℓ, ( ), ∈ [0, 1] 0, otherwise (19)
Multiwavelets in Eq. ( 19) with = 1, 2, 3 are plotted in Fig. [START_REF] Chapelier | Evaluation of a High-Order Discontinuous Galerkin Method for the DNS of Turbulent Flows[END_REF]. At this point, the basis of V 0 , W 0 have been defined. To describe successive subspaces for > 0 we require the mapping ↦ →

2( -) Δ

, with and Δ the center and the size of element , respectively. The basis ℓ , ℓ ∈ V , W are generated by dilation and translation of 0 ℓ0 , 0 ℓ0 ∈ V 0 , W 0 , namely,

ℓ ( ) = 2 Δ ℓ 2( -) Δ , = 0, ..., ; ℓ = 0, ..., ; = 0, ..., 2 -1, (20)
ℓ ( ) = 2 Δ ℓ 2( -) Δ . (21)
Both scaling functions and multiwavelets support extends to the current interval defined by , . That is supp ℓ = supp ℓ = I . Moreover, they are 2 -normalised,

ℓ 2 = ℓ 2 =
1, and share the following orthogonality relations:

ℓ , ℓ I = ℓℓ , ℓ , ℓ I = 0, ℓ , ℓ I = ℓℓ . (22)
Additionally, inherited by Alpert's algorithm, multiwavelets have = ℓ + + 1 vanishing moments, which means that the multiwavelets are orthogonal to polynomials of degree , , ℓ I = 0, ∀ ∈ P -1 . Further details can be found in [START_REF] Hovhannisyan | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws[END_REF].

A strong relation exists between the subspaces defined in the MRA. The multiwavelet subspace W is the orthogonal complement of the the scaling function space V in V +1 [START_REF] Alpert | Adaptive Solution of Partial Differential Equations in Multiwavelet Bases[END_REF], namely,

V +1 = V ⊕ W ; W ⊥ V . ( 23
)
By successive application of Eq. ( 23), a hierarchy of multiwavelet subspaces can be derived:

V +1 = V 0 ⊕ W 0 ⊕ W 1 ⊕ ... ⊕ W -1 . (24)

Coupling of DG and multiwavelets

If we remember the DG method described in Section II, the basis employed for the element-wise solution are built upon Legendre polynomials. As mentioned above, the same occurs with the basis of the scaling function subspace in the multiwavelet formulation. Therefore a direct relation can be established between both approaches. By extending the approximate DG solution in element in Eq. (3) over the entire domain Ω, we have:

ℎ ( , ) = 2 -1 =0 ℓ=0 ℓ ( ) ℓ 2( -) Δ , (25)
where a dyadic mesh composed of = 2 elements has been considered, so that the multiwavelet formulation holds true. It matches the interval previously defined for the multiwavelet formulation. We can now express the DG solution in terms of a combination of scaling functions, namely,

ℎ ( , ) = 2 -1 =0 ℓ=0 ℓ ( ) ℓ = 2 Δ 2 -1 =0 ℓ=0 ℓ ( ) ℓ 2( -) Δ , with Δ = 2 -+1 . (26)
Hence, by comparing Eq. ( 25) and ( 26), the relation between the DG coefficients and the scaling function coefficients is given by

ℓ ( ) = 2 -/2 ℓ ( ). (27)
By using Eq. ( 24), a scaling function subspace can be hierarchically divided into a cascade of multiwavelet subspaces plus a baseline scaling function subspace corresponding to the lowest resolution level. If we apply this to Eq. ( 26), a new representation of the global DG solution is given by:

ℎ ( , ) = ℓ=0 0 ℓ0 0 ℓ0 ( ) + -1 =0 2 -1 =0 ℓ ℓ ( ) , (28)
where

0 ℓ0 = ℎ ( , ), 0 ℓ0 0 0 , ℓ = ℎ ( , ), ℓ ,
are the scaling function and multiwavelet coefficients, respectively. These coefficients can be computed efficiently using the so-called quadrature mirror filter (QMF) coefficients, ℎ (0) ℓ , ℎ (1) ℓ , (0) ℓ , (1) ℓ [START_REF] Geronimo | Alpert Multiwavelets and Legendre-Angelesco Multiple Orthogonal Polynomials[END_REF]. The relations between the coefficients on two consecutive levels are given by:

-1 ℓ = =0 ℎ (0) ℓ ,2 + ℎ (1) ℓ ,2 +1 , (29a)
-1 ℓ = =0 (0) ℓ ,2 + (1) ℓ ,2 +1 . ( 29b
)
Coarser scales of the solution can be obtained by the successive application of Eq. (29a). Furthermore, the multiwavelet details between scales are given by Eq. (29b).

Two-dimensional multiwavelets

In the 2D case the domain in which the scaling functions and multiwavelets operate is Ω = [-1, 1] 2 which is discretised into 2 and 2 elements in the -and -directions, respectively. The results obtained in the previous section are extended in a straightforward manner by using tensor products of the basis. We have the subspace V = span ℓ ( ) ℓ ( ) for the scaling functions, and three components for the multiwavelet subspace:

W ,ℎ = span ℓ ( ) ℓ ( ) (horizontal), W , = span ℓ ( ) ℓ ( ) (vertical), and W , = span ℓ ( ) ℓ ( ) (diagonal)
. Further details can be found in [START_REF] Vuik | Multiwavelet Troubled-Cell Indicator for Discontinuity Detection of Discontinuous Galerkin Schemes[END_REF]. The relation between the DG coefficients and the scaling function coefficients is equivalent to Eq. ( 27) from the previous section, and now is given by

ℓ, = 2 -+ 2 ℓ ; = ( , ), ℓ = (ℓ , ℓ ), = ( , ); ℓ , ℓ = 0, ..., , , = 0, ..., 2 -1 . (30)
By employing the QMF coefficients [START_REF] Geronimo | Alpert Multiwavelets and Legendre-Angelesco Multiple Orthogonal Polynomials[END_REF] the lower-lever scaling and multiwavelet coefficients for each component can be calculated efficiently

-1 ℓ = 1 ¯ , ¯ =0 , =0 ℎ ( ¯ ) ℓ , ℎ ( ¯ ) ℓ , , 2
+ ¯ ; = ( , ), ¯ = ( ¯ , ¯ ) , ( 31a
) ℎ, -1 ℓ = 1 ¯ , ¯ =0 , =0 ( ¯ ) ℓ , ℎ ( ¯ ) ℓ , ,2 + ¯ , (31b)
, -1 ℓ = 1 ¯ , ¯ =0 , =0 ℎ ( ¯ ) ℓ , ( ¯ ) ℓ , ,2 + ¯ , (31c)
, -1 ℓ = 1 ¯ , ¯ =0 , =0 ( ¯ ) ℓ , ( ¯ ) ℓ , ,2 + ¯ . ( 31d
)

IV. h-adaptive algorithm

A multiresolution-based grid adaptivity strategy based on multiwavelets for scalar one-dimensional conservation laws was presented in [START_REF] Hovhannisyan | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws[END_REF]. In that work and in the research that followed [START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF][START_REF] Gerhard | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws: multidimensional case[END_REF][START_REF] Gerhard | Multiwavelet-Based Grid Adaptation with Discontinuous Galerkin Schemes for Shallow Water Equations[END_REF], the problems inherent to a priori and a posteriori error estimates for nonlinear systems of conservation laws were overcome by the use of MRA. That is, the current global solution is simply represented as data on some coarse level plus detail coefficients which embody the fine scale data (see section III). Tresholding then drives the adaptivity process: every multiwavelet coefficient is related to an element in the MRA cascade, nullifying a coefficient is analogous to removing its associated element. Thus, the more details cancelled, the smaller the number of DOFs in the adapted grid. This is the so-called multiwavelet-based grid adaptivity procedure.

In the work presented in [START_REF] Vuik | Multiwavelet Troubled-Cell Indicator for Discontinuity Detection of Discontinuous Galerkin Schemes[END_REF][START_REF] Vuik | Automated Parameters for Troubled-Cell Indicators Using Outlier Detection[END_REF] the same concept of MRA cascading is studied to develop a multiwavelet troubled-cell indicator to identify elements in the vicinity of a shock. However, unlike [START_REF] Hovhannisyan | Adaptive Multiresolution Discontinuous Galerkin Schemes for Conservation Laws[END_REF] in which the entirety of the MRA levels are considered, here the authors base their indicator on the multiwavelet contribution from a single level of the MRA hierarchy.

The relation between the elements in the physical domain and the intervals in the multiwavelet decomposition is bijective. In short, we have one-to-one correspondence. This fact makes the adaptive procedure dependent on the nature of the multiwavelets. To be able to use the dilation and translation properties of the multiwavelets (Eqs. [START_REF] Alpert | Adaptive Solution of Partial Differential Equations in Multiwavelet Bases[END_REF][START_REF] Alhawwary | On the Accuracy and Stability of Various DG Formulations for Diffusion[END_REF]), a sequence of nested dyadic grids is required. Therefore, only meshes which satisfy this condition can be employed. Conversely, in non-Cartesian grids, dilation and translation are no longer available and multiwavelets must be calculated separately for each level and interval of the decomposition. This is an important constraint which increases the cost of the computation. Additionally, due to the duality physical domain/multiwavelet decomposition, difficulties with parallelisation may arise.

As opposed to the above-mentioned MRA strategy, in this work we aim at developing a mesh adaptation algorithm that starts from a coarse solution and proceed with the refinement where required. Traditionally, to deal with this problem, a posteriori error estimators have been used to drive the adaptation, which are computed from the discrete solution [START_REF] Mitchell | A Comparison of hp-Adaptive Strategies for Elliptic Partial Differential Equations[END_REF].

Then we propose a method following the philosophy of an MRA local to the element. More specifically, a MRA is performed independently for every element of the physical domain, in contrast to the global MRA approach which involves using the information from elements on a larger stencil. This approach would allow us to use more general grids (such as locally regular grids) in which elements are not limited to have the same size (but always maintaining a rectiliniar shape). In this case, anisotropic refinement may be applied. Additionally, dealing with each element independently simplifies its use for hp-adaptation. Finally, the excellent parallel properties of the DG method could be perfectly exploited thanks to the local character of the error estimator.

A. Enrichment of DG solution

If we begin our simulations from a coarse grid, the element-wise solution must be enhanced so that the local MRA can be properly applied. This due to the fact that we do not have enough information to start the multiresolution process. We require at least two levels so that we can compare data at two different resolutions and calculate their multiwavelet contribution.

In [START_REF] Dolejší | hp-Discontinuous Galerkin Method Based on Local Higher Order Reconstruction[END_REF], a high order reconstructed solution is measured against the original solution to guide hp-adaptation. The authors verify numerically that the reconstruction better approximates the exact solution compared to the original data. Inspired by this idea, we build the element-wise MRA based on the new locally enriched solution. The new approximation is obtained from the solution within the current element and its neighbours, constructing an overlapping set of patches. It is based on a standard least squares (SLS) reconstruction. A piecewise polynomial is thus constructed from the original solution and extends to the element and its neighbours.

At this point it is important to mention the difference between the meaning we give to parent and children cell. The parent element refers to the current element, denoted by , and the children element refers to the hypothetical isotropic subdivision of the parent element in two (1D) or four (2D) elements, the new elements being denoted by . Only data is stored in a parent/children hierarchy so it can be used in the local MRA. Note that this corresponds to a post-processing step and no actual mesh subdivision actually occurs at this stage.

The reconstruction procedure will be presented first in a one-dimensional setup. The extension to higher dimensions will be described after. The standard least squares (SLS) reconstruction with 2 -norm has been considered in this work. Moreover, analysis of the reconstruction in the parent and children element have been performed as well. In the parent framework we employ a higher order reconstruction on the element by using the solution in the current and neighbouring elements. The reconstruction is then projected on the children with the same degree of the original DG solution. On the other hand, the children framework follows the same philosophy but by performing a direct reconstruction on the children while maintaining the same order of the solution. In total, our study covers two different variations of the SLS. Lastly, an additional method which constructs a new approximation on the children by considering the solution jumps at the interface of the elements is also investigated.

Standard Least-Squares approximation with 2 -norm

The enriched solution is given by a polynomial with degree, = + , where is the degree of the original DG solution and is a parameter that depends on the type (parent or child) of the current element. Therefore, = 0 for a parent element and = 1 when a child element is used. In the latter case, the original discrete domain Ω ℎ is further divided into elements of characteristic size ℎ/2, thus defining Ω ℎ/2 . We define 1 , 2 ∈ Ω ℎ/2 as the left and right half of an element . Consequently, for an element

˜ =          , if element is parent 1 , if element is left child 2 , if element is right child (32)
we consider the sub-domain created by the union of the current and neighboring elements (patch) D = ˜ + ∪ ˜ -. The enriched solution is described by

˜ ℎ ( , ) = ℓ=0 ˜ ℓ ˜ ( ) ℓ ˜ ( ), ∀ ∈ ˜ , ˜ ∈ Ω ℎ , if ˜ = Ω ℎ/2 , if ˜ = 1 , 2 , ∀ > 0 . (33)
The polynomial coefficients ˜ ℓ ˜ are calculated by minimising the error in the least-squares sense

˜ 0 ˜ ( ), . . . , ˜ ˜ ( ) = ∫ D [ ℎ ( , ) -˜ ℎ ( , )] 2 d , (34)
˜ ℓ ˜ = 0, ℓ = 0, . . . . (35)
Solving the optimization problem in Eq. ( 35), we obtain: 00

• • • 0 . . . . . . . . . 0 • • • ˜ 0 ˜ . . . ˜ ℓ ˜ = 0 . . . (36)
where:

= ∫ D ( ) ( ) d , = ∫ D ( ) ℎ ( , ) d , , = 0, . . . , . (37)

Edge reconstruction on children

The reconstruction is performed on the children by considering the jumps at the interfaces between elements. The original discrete domain Ω ℎ is divided into elements of characteristic size ℎ/2, thus defining Ω ℎ/2 . Again, we define 1 , 2 ∈ Ω ℎ/2 as the left and right half of an element . Consequently, for an element

˜ = 1 , if element is left child 2 , if element is right child (38)
we consider the sub-domain resulting by the union of the current element and the shared borders with its neighbours (patch

) D = ˜ ∪ ( ˜ + ∩ ˜ -). The enriched solution is described by ˜ ℎ ( , ) = ℓ=0 ˜ ℓ ˜ ( ) ℓ ˜ ( ), ∀ ∈ ˜ , ˜ ∈ Ω ℎ/2 , ∀ > 0 . ( 39
)
For the left child 1 , the polynomial coefficients ˜ ℓ 1 are obtained by solving the system

0 (-1) • • • (-1) 0 (-1) • • • (-1) 0 ( 0 ) • • • ( 0 ) . . . . . . . . . 0 ( ) • • • ( ) 0 (1) • • • (1) ˜ 0 1 . . . ˜ ℓ 1 = ℎ ( ( 1 -1) + , ) ℎ ( - 1 , ) ℎ ( 0 , ) . . . ℎ ( , ) ℎ ( + 1 , ) (40)
whereas for the right child 2 ,

0 (-1) • • • (-1) 0 ( 0 ) • • • ( 0 ) . . . . . . . . . 0 ( ) • • • ( ) 0 (1) • • • (1) 0 (1) • • • (1) ˜ 0 2 . . . ˜ ℓ 2 = ℎ ( - 2 , ) ℎ ( 0 , ) . . . ℎ ( , ) ℎ ( + 2 , ) ℎ ( ( 2 +1) -, ) (41)
where ℓ refer to the interior Gauss points given by ℓ = 0, ..., and the sub-indices ˜ -, ˜ + denote the values at the left and right edge, respectively.

Two-dimensional reconstruction

In this case we will extend to higher dimensions the standard least squares 2 -norm method given by Eqs. [START_REF] Gartling | A test problem for outflow boundary conditions-flow over a backward-facing step[END_REF][START_REF] Yee | On Spurious Behavior of CFD Simulations[END_REF]. Only the children reconstruction will be presented. With this in mind, parent element will have four children associated. To do that, the original discrete domain Ω ℎ is further divided into elements of characteristic size ℎ/4, thus defining Ω ℎ/4 . Therefore, we have 1 , 2 , 3 , 4 ∈ Ω ℎ/4 and the reconstruction is performed for every child. Consequently, for an element

˜ =              1 , if element is Southwestern child 2 , if element is Northwestern child 3 , if element is Southeastern child 4 , if element is Northeastern child (42)
we consider the sub-domain created by the union of the current and neighboring elements (patch) D = ˜ + ∪ ˜ -. The enriched solution is described by

ũℎ (x, ) = ℓ=1 Ũℓ ˜ ( ) ℓ ˜ (x), ∀ ∈ ˜ , ˜ ∈ Ω ℎ/4 , ∀ > 0 . ( 43
)
We minimise the error in the least-square sense to obtain Ũℓ

( Ũ1 , ..., Ũ ) = ∫ D [u ℎ (x, ) -ũℎ (x, ] 2 d (44) Ũℓ = 0, ℓ = 1, ..., , (45)
which leads to the linear system

Ũℓ = , with = ∫ D ( ) ( ) d , = ∫ D ( ) ℎ ( , ) d , , = 1, ..., . (46)

B. Building local hierarchy via MWs

The reconstruction methods presented in the previous section allow us to construct an enriched solution and establish a local hierarchy between parent (original element) and children (subdivided parent). By means of a MRA within the original element, we will be able to build up a cascade of scales from the enhanced solution in the children up to the parent level. The intensity of the multiwavelet coefficients, which are the building blocks of the local cascading, will be studied in the next section.

To develop a local MRA for every element, a connection between the enriched solution and the multiwavelet analysis is required. In section III the relation between the DG solution and MRA is established by Eq. ( 27) and [START_REF] Taube | hp-Adaptation in Space-Time within an Explicit Discontinuous Galerkin Framework[END_REF] in the oneand two-dimensional settings, respectively. The reconstructed solution is based on Legendre polynomials and behaves very similarly to the DG solution. Therefore the equations still hold true for the reconstructed solution. Consequently, for the one-dimensional analysis we have

˜ ℓ ( ) = 2 -/2 ˜ ℓ ( ) ℓ = 0, ..., ; = 0, ..., 2 -1. ( 47
)
Due to the nature of the enriched solution and the application of the local MRA, a cascade of two levels is attained. Thus we secure the data for the original element (parent) and the next finer level (children). With this in mind, = 1 and Eq. ( 47) is simplified to

˜ 1 ℓ ( ) = 1 √ 2 ˜ ℓ ( ) , ℓ = 0, ..., ; = 0, 1. ( 48
)
where the indices = 0, 1 corresponds to the children 1 , 2 previously defined. The remaining lower-level scaling function and multiwavelet coefficients can be obtained by successively applying the quadrature mirror filter coefficients [START_REF] Geronimo | Alpert Multiwavelets and Legendre-Angelesco Multiple Orthogonal Polynomials[END_REF] developed in Eq. ( 29), which in the particular case of two levels results

˜ 0 ℓ0 = =0 ℎ (0) ℓ ˜ ,0 + ℎ (1) ℓ ˜ ,1 , (49a)
˜ 0 ℓ0 = =0 (0) ℓ ˜ ,0 + (1) ℓ ˜ ,1 . (49b)
We have reached the details, ˜ 0 ℓ0 , with which we opened this section. These will constitute the core of the adaptation process, as will be discussed in the next section. Furthermore, the enriched solution can be now expressed by Eq. ( 28) as a combination of scaling functions and multiwavelets:

˜ ℎ ( , ) = ℓ=0 ˜ 0 ℓ0 0 ℓ0 ( ) + ˜ 0 ℓ0 0 ℓ0 ( ) . (50)
In the two-dimensional setting, Eq. (48) becomes

˜ 1 ℓ, = 1 2 ˜ ℓ ; ℓ = (ℓ , ℓ ), = (
where the indices , = 0, 1 corresponds to the children 1 , 2 , 3 , 4 previously defined. Finally, the QMF coefficients described in Eq. ( 49) should cover the -, -and -components of the multiwavelet coefficients. They read:

˜ 0 ℓ0 = 1 ¯ , ¯ =0 , =0 ℎ ( ¯ ) ℓ , ℎ ( ¯ ) ℓ , ˜ 1 , (0,0)+ ¯ ; = ( , ), ¯ = ( ¯ , ¯ ) , (52a)
˜ ℎ,0 ℓ0 = 1 ¯ , ¯ =0 , =0 (
¯ ) ℓ , ℎ ( ¯ ) ℓ , ˜ 1 , (0,0)+ ¯ , (52b)
˜ ,0 ℓ0 = 1 ¯ , ¯ =0 , =0 ℎ ( ¯ ) ℓ , ( ¯ ) ℓ , ˜ 1 , (0,0)+ ¯ , (52c)
˜ ,0 ℓ0 = 1 ¯ , ¯ =0 , =0 (
¯ ) ℓ , ( ¯ ) ℓ , ˜ , (0,0)+ ¯ . (52d)

C. Thresholding and element marking

By combining the reconstruction techniques with the multiwavelet analysis just discussed, new original error estimators are built. In this section their mathematical description will be developed.

The relation between reconstructed solution and multiwavelets has been established. By imposing an input threshold, the contribution associated to the details derived in Eq. (49b) will be kept or discarded based on their value being higher or lower than the threshold, respectively. Multiwavelet contributions are related to elements. Therefore, discarding a coefficient is equivalent to saying that the element does not require refinement, whereas keeping a coefficient is synonym that further refinement of its associated element is needed. Consequently, the behaviour of the multiwavelet coefficient and the threshold value drive the h-adaptation process. Considering that we have an exact solution to the problem in Eqs.

(1-2), the global 2 -error between the exact solution, ( , ), and the DG approximation, ℎ ( , ), given by solving Eq. ( 4), reads

ℎ = ( , ) -ℎ ( , ) 2 , ∀ ∈ , ∈ Ω ℎ , ∀ > 0 . (53)
Conversely, the global error between the original DG approximation, ℎ ( , ), and the post-enriched DG approximation, ˜ ℎ ( , ), yields

E ℎ = ˜ ℎ ( , ) -ℎ ( , ) 2 , ∀ ∈ , ∈ Ω ℎ , ∀ > 0 . ( 54
)
The input parameter that acts as threshold and controls the adaptation process, , is connected to the error given by Eq. ( 53). The threshold specifies the upper limit error so that we target ℎ to be below this value. Then we would have reached the desired level of accuracy. Therefore, an error estimator which is able to follow the evolution of the error given by the aforementioned equation is a reliable and versatile estimator. Certainly, it does not require third party parameters, its measure is a clear indicator on how the DG approximation behaves, and it is in principle independent on the test case. Our error local estimator is then defined as

= ℓ=0 ˜ 0 ℓ0 0 ℓ0 ( ) 2 = ℓ=0 ˜ 0 ℓ0 2 1/2 , ( 55
)
where the orthonormality relations in Eq. ( 22) simplify the final expression. We now define the global effectivity index as the ratio between the the error given by the estimator, Eq. ( 55), and the error in the approximation, i.e. Eq. ( 53). It reads

eff = ℎ . (56)
where 2 = ∈Ω ℎ 2 . This index is a measure of the quality of the error estimation. The best case scenario is to have an index close to unity. This means that the estimator follows high fidelity the evolution of the error in the approximation. Finally, we describe the set of thresholded elements, that is, elements that are marked for refinement, as Ω ℎ :

Ω ℎ = { ∈ Ω ℎ | > } , (57)
where Ω ℎ is the domain of all elements, as defined in section II, and is the imposed threshold. Furthermore, the change in element size between two neighbouring elements is limited to one. In short, the well known 2:1 condition. This conditions specifies that the difference in size between two elements that share an interface can not be larger than one level of refinement. Due to this constraint, additional elements may be marked for refinement to fulfil that condition. A maximum level of refinement is also specified beforehand to avoid elements of too small size. The simulation is then restarted based on this new mesh. Then, we define the global 2 -error between the exact solution, ( , ), and the DG approximation in the new adapted grid, ℎ ( , ), as follows

ℎ = ( , ) - ℎ ( , ) 2 , with ℎ ( , ) = ℓ=0 ℓ ( ) ℓ ( ), ∀ ∈ , ∈ Ω ℎ , ∀ > 0 . (58)
In the two-dimensional setting, Eq. ( 53) becomes

ℎ = ( , ) -ℎ ( , ) 2 , ∀ ∈ , ∈ Ω ℎ , ∀ > 0 . ( 59
)
and the local error estimator measures the -, -and -components of the multiwavelet contribution in the 2 -norm:

= ℓ ,ℓ =0 ˜ ℎ ℓ ℓ ( ) ℓ ( ) + ℓ ,ℓ =0 ˜ ℓ ℓ ( ) ℓ ( ) + ℓ ,ℓ =0 ˜ ℓ ℓ ( ) ℓ ( ) 2 , (60)
where ˜ ℎ ℓ = ˜ 0,ℎ ℓ0 , ℓ = 0 ℓ0 , and ℓ = 0 ℓ0 for simplicity. Defining

Ψ ℎ ℓ ( , ) = ℓ ( ) ℓ ( ), Ψ ℓ ( , ) = ℓ ( ) ℓ ( ), Ψ ℓ ( , ) = ℓ ( ) ℓ ( ) . (61)
In general, by considering the orthonormality relations in Eq. ( 22), we have

Ψ ℓ , Ψ [-1,1] = , ℓ, ; , = {ℎ, , }, ℓ = (ℓ , ℓ ), = ( , ), ℓ , ℓ , , = 0, ..., . (62)
Consequently, Eq. ( 60) results

=       ℓ ,ℓ =0 ( ˜ ℎ ℓ ) 2 + ℓ ,ℓ =0 ( ˜ ℓ ) 2 + ℓ ,ℓ =0 ( ˜ ℓ ) 2       1/2 . ( 63
)

D. Proposed error estimators

Up to this point, we have recounted two main processes supporting our estimators. In a first step, we have described manifold post-enrichment methods applied to the original DG solution. Subsequent to this step, a local MRA decomposition of the resulting enhanced solution by means of a multiwavelets expansion has been explained. The estimation of the indicators is then calculated by evaluating the 2 -norm of this multiwavelet expansion. The entire procedure has been arranged for one-and two-dimensional settings.

With this framework in mind, we propose three error estimators in total. The difference between them lies in the post-enrichment method employed. The subsequent multiwavelet expansion and its 2 -norm evaluation does not change. Therefore, their terminology is the following:

1) Local multiwavelet parent-based (MW-P) error estimator, defined by MW-P , in which a parent-based 2 -norm SLS reconstruction considering the neighbours is performed. 2) Local multiwavelet child-based (MW-C) error estimator, defined by MW-C , in which a child-based 2 -norm SLS reconstruction considering the neighbours is performed. 3) Local multiwavelet interface-based (MW-I) error estimator, defined by MW-I , in which a child-based reconstruction at the element interface is performed, Eq. (40-41). In this case, parent refers to an element of the current mesh, and children to the post-processing subdivision of the element in two (1D) or four (2D) new elements. Both MW-P and MW-C are grounded on Eqs. (36-37). Lastly, the sub-index denotes that the indicators are local to the element .

Finally, in our study two error estimators extensively tested in [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Discontinuous Galerkin Methods for the Euler and Navier-Stokes Equations[END_REF][START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF] are employed as reference estimators: the small-scale energy density (SSED) [START_REF] Kuru | An Adaptive Variational Multiscale Discontinuous Galerkin Method For Large Eddy Simulation[END_REF] and the spectral decay (SD) [START_REF] Taube | hp-Adaptation in Space-Time within an Explicit Discontinuous Galerkin Framework[END_REF] error estimators. The SSED measures the energy associated with the highest-order modes. It is expressed as:

SSED = ℓ=0 ℓ ( ) ℓ ( ) - -1 ℓ=0 ℓ ( ) ℓ ( ) 2 . ( 64
)
The SD corresponds to the SSED estimator normalised by the total energy within the element. That is:

SD = ℓ=0 ℓ ( ) ℓ ( ) - -1 ℓ=0 ℓ ( ) ℓ ( ) 2 ℓ=0 ℓ ( ) ℓ ( ) 2 .
(65)

V. Numerical results

A. 1D viscous Burgers' equation

A series of steady simulations of the viscous Burgers' equation described in section II are performed. Two classes of initial conditions (ICs) are studied. Firstly, a sharp solution (IC-discontinuity) which simulates the presence of a discontinuity in the middle of the domain is given by ( , 0) =tanh(25 ). The second initial condition is given by a sinusoid and it is representative of a smooth solution (IC-smooth), i.e. ( , 0) =sin(2 ). In the case of a smooth solution, a source term by the method of manufactured solutions is added to the formulation so that the Burgers' equation can reach the sinusoidal solution. To evolve the solution in time from either of the initial conditions up to the steady state, the explicit scheme in Eq. ( 10) is employed. Once the steady state is reached, the error estimators are applied locally to the solution in a post-processing step. is considered in the analysis of each estimator.

Moreover, the global error between the original DG approximation and the post-enriched DG approximation as given by Eq. ( 54) is considered for this test case. There are three versions of this error, corresponding to the three post-enriching methods studied. They are denoted by E P ℎ , E C ℎ , and E I ℎ . Their purpose is to show the difference in behaviour between an estimation based on a multiwavelet expansion, as we propose, and a more traditional approach in which reconstruction and original solution are directly compared, as suggested in [START_REF] Dolejší | hp-Discontinuous Galerkin Method Based on Local Higher Order Reconstruction[END_REF]. We refer to them under the umbrella of post-enriching errors.

Lastly, the small-scale energy density and spectral decay error estimators described in Eqs. ( 64) and (65), respectively, will perform the role of reference indicators. 2a),(2c), and (2e) shows the variation of the exact error in the conservative variable versus the number of degrees of freedom (#DOFs) when uniform ℎ-refinement is performed, ℎ , as well as for the locally ℎ-adapted solution, ℎ , for different degrees of the approximation. In every degree analysed, the error of uniform ℎ-refinement follows the theoretical slope representing the order of the method, as expected. With respect to the adaptive procedure, all the indicators lead to a large decrease in the #DOFs for a given level of accuracy. The SSED and SD show a marginally faster reduction of the #DOFs during the initial refinement steps, but our proposed indicators display a slightly better performance in the last refinement steps, especially for lower degrees of the approximation. For our estimators the savings in #DOFs are around 85% in = 1, Fig. (2a), whereas the reference indicators show a reduction of around 81% in the same degree. The savings are scaled down to 77% and 75%, respectively, when the degree is increased to = 3, Fig. (2e). This shows that the higher the degree, the smaller is the divide between our indicators and the SSED and SD indicators. Indeed, the difference is closed from 4% to 2% when jumping from = 1 to = 3. Overall, the evolution of our error estimators closely resembles the behaviour of the references while showing slightly larger savings at lower degrees. This is an important validation of our proposed estimators.

Finally, adaptive ℎ-refinement driven by the post-enriching errors shows very similar behaviour to the results obtained for our error estimators. This further validates the idea of using reconstruction techniques as an important tool in mesh adaptation. At this point, the question of why using a multiwavelet expansion, which is more computationally expensive, than a direct comparison between the original DG solution and a post-enriched solution may arise. The answer resides in the fact that a multiwavelet expansion gives more information about the solution. This is especially true in higher dimension, where the details are directly given component-wisely. Moreover, studies of the effectivity index (not presented in this work) show a behaviour closer to unity when employing a multiwavelet decomposition. Hence, these reasons motivate the use of a multiwavelet expansion in our error estimators.

We now focus our attention on the ℎ-adapted mesh resulting from the activation of the local MW interface-based, the local MW parent-based, and the local MW child-based estimators. In the previous discussion, the SSED showed a slightly better fit with our indicators compared to the SD, thus in this case we use the SSED as as our only reference indicator. The levels of refinement in the adapted mesh for each estimator are plotted in Fig. (2b), (2d), and (2f) for = 1, = 2, and = 3, respectively. It can be observed how the area surrounding the discontinuity is subjected to a higher level of refinement. This is true for every estimator. When we increase the degree of the solution, the number of refinement levels decreases. This behaviour is expected because we are increasing the spatial resolution by modifying the local polynomial degree and thus fewer elements are required for a prescribed level of accuracy.

The level of refinement of the uniform mesh to fit the prescribed level of accuracy is given by the uniform dashed line. The activation of the MW-P and MW-C estimators leads to an almost identical pattern of refinement centered around the discontinuity, independently of the degree. On the other hand, the MW-I generates a narrower refined region surrounding the discontinuity. For = 3 it even yields a lower refinement level compared to the MW-P and MW-C indicators. In the case of the SSED indicator, the wider refined region below the line speaks of the divide in #DOFs between estimators previously explained. Certainly, it is reduced when the degree is increased but it is always larger compared to our estimators. Again, this is due to the better performance in the last refinement steps of the MW-P, MW-C, and MW-I estimators.

The three novel indicators never surpass the prescribed accuracy line. This behaviour is ideal because it means that their highest refinement level stays below the level of the uniform mesh that imposes the prescribed accuracy. This same behaviour can not be said of the SSED indicator. This indicator exceeds the threshold, obtaining an over-refined adapted mesh. Particularly severe is the behaviour for = 1, surpassing six levels above the threshold. The over-refinement is mitigated in = 3, with only one surplus level. This behaviour is closely related to the effectivity index, as we will discuss in the next paragraph.

A closer look at the evolution of the errors for the adapted and uniform solutions with the initial condition IC-discontinuity and = 3 is illustrated in Fig. [START_REF] Renac | Aghora: A High-Order DG Solver for Turbulent Flow Simulations[END_REF]. In this figure the error of the adapted mesh is compared to the value given by the indicator. If these values show a similar evolution along the refinement process, then the effectivity index, Eq. (56), associated with the indicator is near unity. Values above unity represent > ℎ , whereas values below unity show that < ℎ .

Figure (3a) presents the MW-I indicator. Its effectivity index is between 0.4 < eff < 0.7 and then drops to eff = 0.3 in the last steps of refinement, thus the estimation is always under the approximation error. The same behaviour is observed in Fig. (3b) for the MW-P indicator. However, the first steps of the refinement present a poor effectivity index, which explains the error overshoot in the adapted mesh. Surprisingly, the effectivity index evolves to values close to unity in the last steps of refinement. The MW-C indicator, Fig. (3c), displays a more regular behaviour. Except for the first steps of refinement, the effectivity index remains mostly constant at eff = 2. The estimation evolves along the path of the approximation error, while remaining slightly higher. Last of all, the SSED estimator in Fig. (3d) shows a more erratic behaviour, with a good tracking in the first steps and effectivity ratios progressively deteriorating to eff > 10, which justifies the over-refining in Fig. [START_REF] Chapelier | Inviscid and Viscous Simulations of the Taylor-Green Vortex Flow Using a Modal Discontinuous Galerkin Approach[END_REF].

A further analysis of the effectivity index for the initial conditions IC-smooth and = 3 is performed in Fig. (4). This example shows no difference between uniform and adapted mesh due to the regularity of the solution. In this situation, performing adaptation is not justified because the final adapted grid is nearly indistinguishable from the uniform mesh. However, a few interesting conclusions can be extracted from this analysis. Firstly, the MW-I indicator in Fig. (4a) does not fare too well when the solution is smooth. In this case, the effectivity index progressively drops to eff < 0.1 in the last steps of refinement. Indeed, the influence of the jump of the conservative quantity in the interface between elements for smooth solutions is not well captured by the reconstruction. Thus the disparity between estimation and approximation error. Secondly, the MW-P indicator, Fig. (4b), tracks particularly well the approximation error, with eff = 1 in the mid and last steps of the refinement. However, an unexpected overshoot in the estimator remains above the trajectory of the uniform mesh. This has been likely caused by an under-refinement in the early steps. Observing 3) and ( 4) we conclude that the estimators proposed in this work may constitute a consistent approach for tracking the global error of the DG approximation. They perform substantially better than the SSED estimator, which shows an excessively high value of the effectivity index. In particular, the MW-C estimator displays an uniform behaviour while maintaining an acceptable effectivity index so that it can be reliably used to control the adaptation process. That is why we will be using the MW-C estimator in the study of the two-dimensional backward-facing step flow.

(0, 0) (6ℎ, 0) (14ℎ, 0) (30ℎ, 0) 5ℎ 60ℎ ℎ Mesh A Mesh B Mesh C Mesh #Elts #DOFs #Elts- #DOFs- = 1 = 2 = 3 = 1 = 2 = 3 A 2144

B. Laminar flow over a backward-facing step

In this section the numerical solution of a two-dimensional steady laminar incompressible flow over a backward-facing step is employed as starting point to prove the validity of our novel error estimator to perform mesh adaptation.

The geometry of the domain is shown at the top of Fig. [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF]. The inlet boundary is located five step heights upstream of the step and the outlet or exit boundary is placed at 60 step heights downstream of the step. The height of the inlet channel is equal to the dimension of the step, and the channel height in the expanded region (downstream of the step),

, is twice the height of the step, ℎ. Therefore, the expansion ratio of the backward-facing step results in /ℎ = 2. The Reynolds number of the problem is set at Re = 800. The Reynolds is defined as Re = where is the inlet mean velocity, i.e. two thirds of the maximum inlet velocity.

With respect to the boundary conditions, we impose at the inlet boundary a fully developed plane Poiseuille flow so that the velocity profile at the entrance of the domain is parabolic. An extra inlet channel has been considered to negate the influence of the step in the upstream flow region [START_REF] Barton | The Entrance Effect of Laminar Flow Over a Backward-Facing Step Geometry[END_REF]. The study [START_REF] Erturk | Numerical Solutions of 2-D Steady Incompressible Flow Over a Backward-Facing Step, Part I: High Reynolds Number Solutions[END_REF] concludes that the length of the inlet channel must be at least five times the height of the step. Thus we have sized our inlet channel accordingly. At the exit boundary, a non-reflecting boundary condition is imposed such that the velocity profile of the numerical solution at the exit boundary matches the analytical parabolic profile of a Poiseuille flow. Indeed to that end, the exit boundary is placed at a distance sufficiently away from the step so that the flow becomes fully developed. The studies [START_REF] Keskar | Computations of a Laminar Backward-Facing Step Flow at Re=800 with a Spectral Domain Decomposition Method[END_REF][START_REF] Gartling | A test problem for outflow boundary conditions-flow over a backward-facing step[END_REF], where a backward-facing step flow at Re = 800 is also investigated, place the exit boundary at 60 step heights downstream from the step. Therefore, we have sized our expanded channel as such. Finally, wall boundary conditions are imposed on the upper and lower surface of the channel, as well as on the surface of the step. Further details are reported in Fig. [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF].

To mesh the domain we have employed four different Cartesian uniform grids, depending on the order of the numerical solution at hand. At the bottom row of Fig. [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF], the grids are described. An overview of the first three meshes, from coarser to finer, is displayed as well. The grids are named A to D, and ordered by increasing resolution. The color nomenclature puts in the same group those grids with equivalent resolution with respect to the number of degrees of freedom (#DOFs). That means that their numerical solution should be comparable. For those very coarse grids (mesh A and B), the area following the inlet has been locally refined so that we are able to properly capture the parabolic inflow. The study [START_REF] Yee | On Spurious Behavior of CFD Simulations[END_REF] shows that when a coarse grid is employed in the backward-facing step flow at high Reynolds numbers, a spurious oscillating numerical solution is obtained and the steady state can not be reached. The work [START_REF] Erturk | Numerical Solutions of 2-D Steady Incompressible Flow Over a Backward-Facing Step, Part I: High Reynolds Number Solutions[END_REF], in which an interval of Re = [100, 3000] is investigated by solving the flow using a second order finite difference method, employs a very fine mesh so that convergence to steady state can be achieved. We have used a DG method in our study and the convergence problems reported in references [START_REF] Erturk | Numerical Solutions of 2-D Steady Incompressible Flow Over a Backward-Facing Step, Part I: High Reynolds Number Solutions[END_REF][START_REF] Yee | On Spurious Behavior of CFD Simulations[END_REF] were not encountered, even for solutions with low polynomial degrees.

To evaluate the quality of our numerical solution, we have selected three streamwise locations along the expanded channel, Fig. [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF]. They include the lower and upper recirculation bubbles, and an overview of the developed flow far away downstream. The idea is to extract the profiles of the relevant physical quantities along the vertical direction of the main channel. We will consider the profiles of the horizontal and vertical components of the velocity, given by and . The profile of the vorticity, defined as = / -/ , is also included in the analysis. These profiles are examined and compared to the results from [START_REF] Erturk | Numerical Solutions of 2-D Steady Incompressible Flow Over a Backward-Facing Step, Part I: High Reynolds Number Solutions[END_REF] at the same streamwise locations.

In this study a grid of 101 uniform elements along the vertical direction of the expanded channel is used. Their scheme is second order accurate. Thus we count 202 DOFs along the -direction. To perform a fair comparison, we have computed a second order numerical solution with 256 DOFs along the same direction. However, we have set side by side this detailed outcome with coarser solutions of our own (not reported in this study) and found very little difference in the profiles analysed. This means that we reached mesh convergence, and the extra DOFs do not significantly change the solution. Therefore, to rely on the grids reported in Fig. ( 5), we will be using a numerical solution with = 3 and computed on mesh C. In this manner, we report 128 DOFs along the -direction of the main channel.

The results are shown in Fig. [START_REF] Harten | Multiresolution Representation of Data: A General Framework[END_REF]. We observe that our computed profiles agree well with those of the literature for every physical quantity analysed, even though we use nearly half the #DOFs. There are a small discrepancy in the maximum values of the -and -velocities. We believe this is due to the fact that in our simulations the maximum -velocity of the parabolic profile is slightly higher right at the fall of the step than at the upstream inlet. In [START_REF] Erturk | Numerical Solutions of 2-D Steady Incompressible Flow Over a Backward-Facing Step, Part I: High Reynolds Number Solutions[END_REF] it is shown that for Re = 800, these two profiles must be nearly identical. We believe that by using a compressible solver (for an incompressible flow) we would need to further adjust the BC to mitigate the small increase in the horizontal velocity.

After validating our numerical results on uniform grids, the idea is to start from a coarse mesh (such as mesh A or B) so that our error estimator may detect the areas in which the solution is underresolved.

Once the steady-state solution is reached, a post-processing step is carried out in which an enhanced solution is constructed from the converged numerical solution in a local manner, by following the two-dimensional 2 -norm standard least squares method proposed in section IV. In this case, the least squares method is constructed from the contribution of the current cell and neighboring children. This process has already been described in Eqs. (44 to 45). The multiwavelet transform is then applied to the new enriched solution, decomposing the solution into a coarse and a fine resolution components in the horizontal, vertical and diagonal directions. This process may be performed for every conservative variable u = ( , v, ) or any other derived quantity, such as the vorticity. For this study we select the horizontal and vertical component of the momentum vector, which for this configuration constitutes a relevant quantity representative of the behaviour of the overall solution.

The main goal of this work is to illustrate the performance of our error estimator in the framework of isotropic ℎ-adaptation. Therefore, the contribution of each direction is merged into one single quantity that will be assigned to the analysed element. This is done by applying Eq. ( 63) in an element-wise manner. On the other hand, if anisotropic ℎ-adaptation were to be performed, each component would be studied separately. 5), we have selected grids with similar #DOFs along the -direction to obtain equivalent results when the degree of the numerical solution is changed. That is, mesh A for = 2, 3 and mesh B for = 1.

We observe from Fig. (7a) that for = 1 the estimation reports high values at the inlet channel and the upper region immediately after the step. This means that the velocity profile is not well resolved upstream the step and a finer grid is required. On the other hand, the recirculation bubble caused by the separation of the flow at the step corner displays low estimator levels, which can be explained by the regularity of the solution in that region.

If we move our attention to = 2, 3, we can see that the velocity profile at the entrance of the channel is somewhat better captured. This is reflected by a decrease of the indicator value in that region. Somehow moderate estimation values with = 2 (when a parabola should be exactly captured) are likely due to the small influence of the step upstream and the fact that the numerical solution is underresolved on this mesh. Then the parabolic profile might be not fully captured by the numerical solution. The recirculation region continues to display low indicator values as well. However, -10 0 10 vorticity,

(h) = 2.
-10 0 10 vorticity, However, they fully disappear when the order of the numerical solution is increased. Adaptation in the highest order solution, Figs. (8c); (8f); and (8i), reports excellent agreement with the reference.

(i) = 3.
We have inspected the profiles of several physical variables at the step. This analysis has shown very promising results when adaptation is activated. Moreover, we can conclude that the static pressure is the most sensitive quantity to adaptation. Therefore, we will use this variable to evaluate how adaptation influences the behaviour of the flow downstream.

In Fig. [START_REF] Roussel | A Conservative Fully Adaptive Multiresolution Algorithm for Parabolic PDEs[END_REF] we plot the pressure profiles for = 1, 2, 3 at the previously defined downstream locations of = 6ℎ, = 14ℎ, and = 30ℎ. In ascending order, each row refers to one location.

The combination of coarse or medium grids (dotted and black solid line, respectively) and low order solutions causes the pressure to oscillate compared to the reference (marker). This is clearly observed at = 6ℎ, Fig. . This time only due to underresolving, as we are away from the singularity. When we increase the order of the solution the oscillations are lessened.

The adapted pressure (solid red line) reports as good behaviour as the medium grid (with slight oscillations) and becomes closer to the reference solution at the location nearest to the step, Figs. (9a); (9b); (9c). At farther locations the adaptations fares slightly worse than the medium grid for low order solutions, Figs. (9d) and (9g). However this trend is reversed with ever increasing order. The highest order adapted solutions in Figs. (9c), (9f), and (9i) show a very good match with the reference, slightly improving coarse and medium grids. Interestingly, we suspect that by refining around the upstream singularity, the adapted solution downstream gets somewhat closer to the reference solution. Thus correctly capturing the singularity is especially relevant.

We conclude our study with the analysis of the computational times for the different grids evaluated. The survey presented in Table (2) records the times for every numerical simulation of the backward-facing step flow performed in the study. Simulations have been performed on 96 cores distributed on four nodes of 24 cores each by combining hybrid MPI/OpenMP strategies. If we return to the data in Table [START_REF] Chapelier | Evaluation of a High-Order Discontinuous Galerkin Method for the DNS of Turbulent Flows[END_REF], we obtained #DOFs savings above 90% when performing adaptation. Speedup times for the final adapted grids vary depending on the order of the adapted solution. We reached almost 40 times faster solutions with = 3, nearly 20 times with = 2, and 80 times with = 1. Estimation, marking and refining algorithm times have also been considered when performing adaptation. It is worth mentioning that they never constitute more than 5% of the total computational time. Computational times need to be interpreted with caution because they are hard to measure consistently and are subjected to many variables not always fully understood. However, from this data we can conclude that by activating adaptation we obtain an important reduction both in the #DOFs and simulation times.

VI. Conclusions

In this work we have presented a novel error estimator based on an local multiwavelet analysis of a DG solution which is subjected to an post-enrichment process. This operation is performed so that the multiwavelet expansion is able to extract meaningful information from the improved approximation of the original DG solution. The new information is represented in the form of multiwavelet coefficients, on which is based the computation of the error estimator. Simulations of the one-dimensional viscous Burgers equation show promising results. An important decrease in the number of degrees of freedom (#DOFs) is observed when the adaptive process is activated. A comparison with the modal error estimators developed in [START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Discontinuous Galerkin Methods for the Euler and Navier-Stokes Equations[END_REF][START_REF] Naddei | A Comparison of Refinement Indicators for -adaptive Simulations of Steady and Unsteady Flows Using Discontinuous Galerkin Methods[END_REF] shows that, for our proposed estimator, fewer elements are adapted for a similar global level of the error. Interestingly, the effectivity index of the proposed method is close to unity. Supported by the encouraging results of the one-dimensional simulations, we extend our estimator to higher dimensions. The early tests are performed on a steady laminar backward-facing step flow at Re = 800. By using a coarse initial mesh we observe high values of the estimator in regions where the flow is underresolved or a singularity manifests. In particular, the areas along the separated shear layer and around the corner of the step, where a sharp velocity gradient exists due to presence of the geometrical jump, are marked for refinement. We can therefore conclude that the relevant regions are successfully detected and the grid adapted accordingly. By performing adaptation, we drastically reduce the #DOFs and computational times while maintaining a analogous behaviour and similar accuracy with respect to the reference grid. Particularly, adaptation is well suited to deal with singularities, where the gains are most noticeable. Most importantly, for a similar #DOFs, adaptation seems to perform better the higher the order we use. This result highlights the importance of high-order methods in adaptation. We are currently developing the extension to isotropic hp-adaptation by monitoring the rate of decay of the multiwavelets' vanishing moments. Future research is also planned on the analysis of the multiwavelet components in the -, -and -directions, which will open the door to the development of an anisotropic adaptation algorithm.

Fig. 1

 1 Fig. 1 The scaling functions (left, up to = 2) are given by Legendre polynomials in [-1, 1]. Different plots of the multiwavelets are displayed for = 0 (2nd from left), = 1 (2nd from right), and = 2 (right). The building algorithm is described in [19].

  , ); ℓ , ℓ = 0, ..., , , = 0, ..., 1 ,

Fig. 2 1D

 2 Fig. 2 1D viscous Burgers' equation with IC-discontinuity. Comparison of different error estimators. The terminology of the estimators reads 1 = MW-I , 2 = MW-C , 3 = MW-P , 4 = SSED , and 5 = SD .

Fig. 3

 3 Fig. 3 1D viscous Burgers' equation with IC-discontinuity, = 3. Different estimators compared. The approximation error is plotted for the uniform and adapted mesh . The global estimation is given by .

Fig. 4 1D

 4 Fig. 4 1D viscous Burgers' equation with IC-smooth, = 3. Different estimators compared. The approximation error is plotted for the uniform and adapted mesh . The global estimation is given by .

Fig. 5

 5 Fig. 5 Backward-facing step domain and boundary conditions (top) and examples of the uniform grids employed (bottom left). The Table (bottom right) accounts the number of elements and degrees of freedom, both in total and along the -direction of the expanded channel.

Fig. 6

 6 Fig. 6 Profiles for different physical quantities at three different streamwise locations along the expanded channel. The profiles have been extracted from our reference numerical solution and compared with the literature at the same locations.

( a )

 a Value of the MW-C estimator for = 1, 2, 3; in logarithmic scale. (b) Final mesh adaptation for = 1, 2, 3 around the region of the step.

Fig. 7

 7 Fig. 7 Overview of error estimation and mesh adaptation.

  Fig. (7a) shows a close-up view of the distribution of the estimation of the error in the region around the step. The value of the estimator is highly dependent on the degree of the approximation and the number of elements. As described by the color nomenclature in Fig. (

Fig. 8 Fig. 9

 89 Fig. 8 Profiles for different physical quantities at step, = 0, where a singularity manifests. Uniform (capital letters) and ℎ-adapted meshes (capital plus Roman numerals) are compared.

  (9a), and = 14ℎ, Fig. (9g)

Table 1 #DOFs savings between uniform and adapted grids.

 1

	Mesh		#DOFs		Mesh #DOFs	Mesh	#DOFs			#DOFs
		= 1	= 2	= 3	B	= 1	A	= 2	= 3	= 1	= 2	= 3
	A
	B		72 864	129 536			I	43 785 35 888		39.9% 72.3%
	C	128 000	288 000	512 000	I	89 000	II	51 075 36 512	30.5% 82.3% 92.9%
	D	512 000	1 152 000		II	101 456	III	51 183		80.2% 95.6%
	E	2 048 000			III	101 708				95.1%

Table 2 Simulation time speedups between uniform and adapted grids.

 2

	Mesh		Time (s)		Mesh Time (s)	Mesh	Time (s)	Time speedups
		= 1	= 2	= 3	B	= 1	A	= 2	= 3	= 1	= 2	= 3
	A
	B		319	463			I	727	230		0.4	2.0
	C	1240	1520	5817	I	1581	II	1218	154	0.8	1.2	37.8
	D	11 248 23 061		II	1038	III	1420		10.8	16.2
	E	81 398			III	1016				80.1

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 675008 for the years 2017 to 2019 and from ONERA for 2020 onwards.

the most interesting aspect is the concentration of elevated estimator values following the separated shear layer and in the vicinity of the step corner. Indeed, a strong velocity gradient appears due to presence of the geometrical jump.

The estimation shown in Fig. (7a) sheds light on the local behaviour of the numerical solution and allows us to detect these regions in which the flow may require a higher order in the approximation or a smaller size of the element. The latter is the main focus of this work. Consequently, the estimator will flag those elements which exceed the value of threshold, . Every marked element is then isotropically divided into four new elements and a new grid is generated. The previously converged numerical solution obtained in the preceding coarser grid is then projected on the new adapted grid and the discrete problem is solved again until the steady state is reached. Once the new numerical solution is achieved, the error estimation process is repeated on the new approximation to update the value of the indicator on the adapted mesh. This recursive procedure continues until the element-wise estimator is below the given threshold.

In Fig. (7b) the final adapted grids resulting from our adaptive algorithm are shown for different degrees of the numerical approximation.

The adapted grid associated to the lowest degree = 1 exhibits the largest number of new elements overall compared to the initial mesh B. These new elements are mostly clustered along the inlet channel and along the upper half of the expanded channel. The separated shear layer is also further refined. In comparison, the recirculation region remains virtually unrefined. These observations are consistent with the error indicator maps shown in Fig. (7a).

On the other hand, the final adapted grid associated with = 3 remains largely unaltered in comparison with the initial mesh A. Only the shear layer and the area around the corner are refined. Then again, the results are in line with the estimated error distribution observed in Fig. (7a). Regarding the final adapted grid obtained for the quadratic approximation = 2, it appears to be in between those obtained for = 1 and = 3. Certainly, the upper part of the domain has been uniformly refined. Though the refinement is coarser than for = 1. As expected, further refinement is performed along the shear layer and the region surrounding the corner.

A full analysis of the grids involved in the adaptation process is presented in Table [START_REF] Chapelier | Evaluation of a High-Order Discontinuous Galerkin Method for the DNS of Turbulent Flows[END_REF]. Only uniform grids with element size at least equal to the minimum element size of the adapted grid (symbolised by Roman numerals) are displayed. One additional uniform grid not reported before (mesh E) is added to the list for = 1. Moreover, we extend = 2 for mesh D as well. The adapted grids are separated by the order of the solution. That is, for = 1 we go through three levels of refinement to eventually obtain mesh B(III). Three refinement steps are also performed for = 2, reaching mesh A(III). Finally, with = 3 we achieve two levels of refinement, obtaining mesh A(II). The idea is to compare the #DOFs of the uniform grids against those of the adapted grids with the same smallest element size. With this in mind, Table [START_REF] Chapelier | Evaluation of a High-Order Discontinuous Galerkin Method for the DNS of Turbulent Flows[END_REF] reports an important reduction in the #DOFs for every successive adapted grid in comparison with their uniform counterparts. The longer the adaptation goes on, the larger the savings. Moreover, we report a faster rate the higher the order of the solution. The savings for the final adapted grids are above 90%.

Certainly the adaptation process has considerably reduced the #DOFs required to solve in the computations. Now we need to verify that the behaviour and accuracy of the numerical solution is not deteriorated despite the downscaling of the #DOFs.

To do so, we investigate how the actual physical flow is affected by the adaptation, and if the adapted solution is as accurate as the one computed on the uniform grids.

Firstly, we evaluate the area around the step, where the singularity occurs due to the geometrical jump. Profiles of the static pressure, -velocity, and vorticity are reported in Fig. [START_REF] Müller | Multiresolution Schemes for Conservation laws[END_REF]. Each row corresponds to one of these quantities while columns refer to the order of the numerical solution. The uniform reference grid (mesh D or C, depending on the order) is denoted by the marker. They constitute the fine grids to match. Coarse and medium uniform grids (mesh A, B or C; depending on the order) are also considered to monitor the effect of the adaptation.

The most dramatic oscillations are reported in the coarser grids (dotted line). This wavering behaviour is especially pronounced in the low order solution, particularly for the pressure, Fig. (8a), and vorticity, Fig. (8g). This is due to two factors. First, the solution is underresolved and two, the presence of the singularity. The use of a medium grid (black solid line) somehow removes the problem of underresolving. Certainly, pressure stops oscillating outside the singularity and the -velocity is closer to reference, Fig. (8d). However, the vorticity still reports fluctuations.

If we increase the order of the solution the underresolving is mitigated, but the effect of the singularity remains for both coarse and medium grids, however with smaller amplitude. Pressure, Fig.