¸ois Vilar
email: francois.vilar@umontpellier.fr

R Émi

A POSTERIORI LOCAL SUBCELL CORRECTION OF HIGH-ORDER DISCONTINUOUS GALERKIN SCHEME FOR CONSERVATION LAWS ON TWO-DIMENSIONAL UNSTRUCTURED GRIDS

Keywords: a posteriori correction, subcell correction, arbitrary high-order, DG subcell FV formulation, positivitypreserving scheme, hyperbolic conservation laws, subcell conservative scheme

. The technique is based on the reformulation of DG scheme as a finite volume (FV) like method through the definition of some specific numerical fluxes referred to as reconstructed fluxes. High-order DG numerical scheme combined with this new local subcell correction technique ensures positivity preservation of the solution, along with a low oscillatory and sharp shocks representation.

The main idea of this correction procedure is to retain as much as possible the high accuracy and the very precise subcell resolution of DG schemes, while ensuring the robustness and stability of the numerical method, as well as preserving physical admissibility of the solution. Consequently, an a posteriori correction will only be applied locally at the subcell scale where it is needed, but still ensuring the scheme conservation. Practically, at each time step, we compute a DG candidate solution and check if this solution is admissible (for instance positive, non-oscillating, . . . ). If it is the case, we go further in time. Otherwise, we return to the previous time step and correct locally, at the subcell scale, the numerical solution. To this end, each cell is subdivided into subcells. Then, if the solution is locally detected as bad, we substitute the DG reconstructed flux on the subcell boundaries by a robust first-order numerical flux. And for subcell detected as admissible, we keep the high-order DG reconstructed flux which allows us to retain the very high accurate resolution and conservation of the DG scheme. As a consequence, only the solution inside troubled subcells and its first neighbors will have to be recomputed, elsewhere the solution remains unchanged. Another technique blending in a convex combination fashion DG reconstructed fluxes and first-order FV fluxes for admissible subcells in the vicinity of troubled areas will also be presented and prove to improve results in comparison to the original algorithm introduced in [52]. Numerical results on various type problems and test cases will be presented to assess the very good performance of the design correction algorithm.

1. Introduction. This paper is devoted to the extension of the a posteriori local subcell correction (APLSC) method introduced in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] to the two-dimensional unstructured case. It is well known that the Cauchy problem for any non linear hyperbolic problem leads to handling, in the generic situation, solutions with discontinuities. In numerical simulations, this aspect has to be dealt with, and this issue has been one of the main questions to address for designing numerical methods for hyperbolic problems. Another core problem is the one of accuracy. The last one is that the problem is well set only if the solution belongs to an admissibility set called invariant domain. For example, for the gas dynamics Euler equations, the density and the internal energy must stay positive. This is particularly difficult to address, simultaneously, these three questions together, because those constraints are antagonistic. These questions has been at the center of algorithmic developments since decades, one may mention [START_REF] Boris | Flux-corrected transport. I: SHASTA, a fluid transport algorithm that works[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. V-A second-order sequel to Godunov's method[END_REF][START_REF] Harten | Uniformly high order essentially non-oscillatory schemes[END_REF][START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments[END_REF] and the reference herein.

The Discontinuous Galerkin (DG) method is one of the most widely used numerical scheme, especially in the context of computational fluid dynamics. Cockburn, Shu et al in a series papers [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems[END_REF][START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case[END_REF][START_REF] Cockburn | The Runge-Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems[END_REF] have paved the road to efficient methods for fluid dynamics. DG methods allow to reach any arbitrary order of accuracy, while keeping the stencil compact, along with other good properties such as built-in entropy stability and hp-adaptivity. However, though recent progress have been done [START_REF] Shu | TVB uniformly high-order schemes for conservations laws[END_REF], designing method that are oscillation free and compliant with the invariant domain is not a trivial matter. Controlling spurious oscillations has been studied in many papers, among which [START_REF] Biswas | Parallel adaptive finite element methods for conservation laws[END_REF][START_REF] Burbeau | A problem-independent limiter for high-order Runge Kutta discontinuous Galerkin methods[END_REF][START_REF] Krivodonova | Limiters for high-order discontinuous Galerkin methods[END_REF][START_REF] Yang | A parameter-free generalized moment limiter for high-order methods on unstrucured grids[END_REF][START_REF] Park | Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids[END_REF][START_REF] Kuzmin | Slope limiting for discontinuous Galerkin approximations with a possibly non-orthogonal Taylor basis[END_REF][START_REF] Li | A new vertex-based limiting approach for nodal discontinuous Galerkin methods on arbitrary unstructured meshes[END_REF]. Staying in the invariant domain has also been made possible, see for example [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments[END_REF][START_REF] Wu | Geometric quasilinearization framework for analysis and design of bound-preserving schemes[END_REF][START_REF] Guermond | Entropy viscosity method for high-order approximations of conservation laws[END_REF]. However, this is often achieved to the cost of enlarging the width of the discontinuous patterns, and some time a lost of accuracy.

Other methods have recently gained in popularity, the so-called subcell techniques. Here the idea is to subdivide the bad cells, and to adopt a special procedure with the hope of curing the negative aspects of the original scheme. Examples of this strategy can be found in [START_REF] Huerta | A simple shock-capturing technique for high-order discontinuous Galerkin methods[END_REF][START_REF] Casoni | One-dimensional shock-capturing for high-order discontinuous Galerkin methods[END_REF][START_REF] Sonntag | Shock capturing for discontinuous Galerkin methods using finite volume subcells[END_REF][START_REF] De La Rosa | Hybrid DG/FV schemes for magnetohydrodynamics and relativistichydrodynamics[END_REF]. For example, [START_REF] Huerta | A simple shock-capturing technique for high-order discontinuous Galerkin methods[END_REF], the authors use a convex combination between high-order DG schemes and first-order Finite-Volumes (FV) on a subgrid, allowing them to retain the very high accurate resolution of DG in smooth areas and ensuring the scheme robustness in the presence of shocks. Similarly, in [START_REF] Sonntag | Shock capturing for discontinuous Galerkin methods using finite volume subcells[END_REF][START_REF] De La Rosa | Hybrid DG/FV schemes for magnetohydrodynamics and relativistichydrodynamics[END_REF], after having detected the troubled zones, cells are then subdivided into subcells, and a robust first-order finite volume scheme is performed on the subgrid in troubled cells. Alternatively, some robust high-order scheme as MUSCL or WENO could either be used to avoid too much accuracy discrepancy. Note that in general this methods are tuned for Cartesian meshes.

Another approach which is worth to be mentioned is the so-called MOOD technique, [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF][START_REF] Diot | Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The MOOD method in the three-dimensional case: very-high-order finite volume method for hyperbolic systems[END_REF]. Through this procedure, the order of approximation of the numerical scheme is locally reduced in an a posteriori sequence until the solution becomes admissible, i.e oscillation free and the solution lives in the invariant domain. In [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws[END_REF][START_REF] Dumbser | A simple robust and accurate a posteriorisub-cell finite volume limiter for the discontinu-ousGalerkin method on unstructured meshes[END_REF], a subcell FV technique similar to the one presented in [START_REF] Sonntag | Shock capturing for discontinuous Galerkin methods using finite volume subcells[END_REF] has been applied to the a posteriori paradigm. In practice, if the numerical solution in a cell is detected as bad, the cell is then subdivided into subcells and a first-order finite volume, or alternatively other robust scheme (second-order TVD FV scheme, WENO scheme, . . . ), is applied on each subcell. Then, through these new subcell mean values, a high-order polynomial is reconstructed on the primal cell. This correction procedure has the benefit to be very simple and robust, and is able to preserve the high accuracy of DG schemes in smooth areas.

In all these aforementioned limitation techniques, a priori and a posteriori, the high-order DG polynomial is either globally modified in the cell, or even discard as it is in the (H)WENO limiter or any a posteriori correction technique in the troubled cells. Since one of the main advantage of high-order scheme is to be able to use coarse grids while still being very precise, one can see that there is a waste of information here, as well as unnecessary computational effort made. This is particularly the case in the vicinity of discontinuities since the polynomials are globally modified. This problem was addressed in the one-dimensional case in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. This new technique relies on the reformulation of DG schemes as a FV-like scheme defined on a subgrid. First, as the number of subcells matches the dimension of the polynomial solution space, the numerical solution inside a cell can be uniquely defined by either a high degree polynomial, or as a piecewise constant solution through its different subcell mean values; the connection between the two being done via a projector. Second, by means of particular basis functions introduced in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], which are nothing but the L 2 projection over the polynomial space of the subcells indicator function, the DG volume and boundary terms can be rewritten as flux differences. Such theoretical analysis is relatively simple in one dimension in space, but much more challenging in two dimensions.

The format of the paper is as follow. Extending the theoretical analysis introduced in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] and using ideas from [START_REF] Abgrall | Some Remarks About Conservation for Residual Distribution Schemes[END_REF], we first reinterpret unstructured grid DG scheme as a subgrid FV-like scheme, through the definition of particular fluxes that we referred to as reconstructed fluxes. These DG reconstructed fluxes are analytically computed, and the analysis shows how those fluxes are connected to the interior polynomial flux and the jump at the cell interface between the interior flux and the DG numerical flux. Let us emphasize that reformulation of DG scheme as a subcell finite-volume method can be performed regardless the form of the element: this is thus not limited to quadrilateral nor triangular cells, but can be done on general polygonal elements. This theoretical part is done in section 2, where a discussion on the type of subcells is also provided. Using this equivalent formulation, we can proceed by means of a MOOD paradigm as follow: at each time step, we compute a DG candidate solution and check if this solution is admissible. If it is the case, we go further in time. Otherwise, we return to the previous time step and correct locally, at the subcell scale, the numerical solution. In the subcells where the solution was detected as bad, we substitute the DG reconstructed flux on the subcell boundaries by a robust first-order numerical flux. And for subcells detected as admissible, we keep the high-order reconstructed flux which allows us to retain the very high accurate resolution and conservation of DG schemes. Consequently, only the solution inside troubled subcells and their first neighbors will have to be recomputed. Elsewhere, the solution remains unchanged. This correction procedure is then extremely local. Another technique blending in a convex combination fashion DG reconstructed fluxes and first-order FV fluxes for admissible subcells in the vicinity of troubled areas will also be presented and prove to improve results in comparison to the original algorithm. This procedure and all the technical details are described in section 3. Section 4 provides numerical results, and a conclusion follows.

2. Discontinuous Galerkin method as a subcell Finite Volume scheme. This section is devoted to the demonstration of the equivalency between DG schemes and a FV-like method on a subgrid provided the definition of particular fluxes. To remain as simple as possible, two-dimensional Scalar Conservation Laws (SCL) will be considered. Let u = u(x, t), for x ∈ ω ⊂ R 2 and t ∈ [0, T ], be the solution of the following system, (2.1a)

(2.1b)    ∂ t u(x, t) + ∇ x F (u(x, t)) = 0, (x, t) ∈ ω × [0, T ], u(x, 0) = u 0 (x), x ∈ ω,
where u 0 is the initial data and F (u) the 2D flux function. For the subsequent discretization, let us introduce the following notation. Let {ω c } c be a generic partition of the domain ω into non-overlapping cells, with |ω c | being the size of ω c . We also partition the time domain in intermediate times (t n ) n with ∆t n = t n+1 -t n the n th time step. In the DG frame, the numerical solution is considered piecewise polynomial over the domain, and hence developed on each cell onto P k (ω c ), the set of polynomials of degree up to k defined on cell ω c . This space approximation theoretically leads to a (k + 1) th space order accurate scheme. Let u c h be the restriction of u h , the piecewise polynomial approximation of the solution u, over the cell ω c

u c h (x, t) = N k m=1 u c m (t) σ c m (x), (2.2)
where the u c m are the N k successive components of u c h over the polynomial basis {σ c m } m . We recall that in the two-dimensional case

N k = (k+1)(k+2) 2 . The coefficients u c m present in (2.
2) are the solution moments to be computed through a local variational formulation on ω c . To this end, one has to multiply equation (2.1a) by ψ ∈ P k (ω c ), a polynomial test function, and integrate it on ω c . By means of an integration by parts and substituting the solution u by its approximated polynomial counterpart u c h , one gets

ωc ∂ u c h ∂t ψ dV = ωc F (u c h ) ∇ x ψ dV - ∂ωc ψ F n dS, ∀ ψ ∈ P k (ω c ). (2.3)
The numerical solution u c h is then the unique polynomial function defined in P k (ω c ) satisfying equation (2.3) for all function ψ ∈ P k (ω c ). In (2.3), the numerical flux function F n , in addition to ensure the scheme conservation, is the cornerstone of any finite volume or DG scheme regarding fundamental considerations as stability, positivity and entropy among others. In this context, this numerical flux is defined as a function of the two states on the left and right of each interface

F n = F u -, u + , n , (2.4) with u -= lim ϵ → 0 + u c h (x i -ϵ n, t) and u + = lim ϵ → 0 + u v h (x i +ϵ n, t)
, where ω v is a face neighboring cell of ω c , while x i and n respectively stand for a point and the unit outward normal of the separating interface. From now on, we refer V c to as the set containing the face neighboring cells of ω c . This function is generally obtained through the resolution of an exact or approximated Riemann problem. In the remainder of this paper, for sake of simplicity, we make use of the very well-known global Lax-Friedrichs numerical flux which reads

F(u, v, n) = F (u) + F (v) 2 n - γ 2 (v -u), (2.5)
where γ = sup w (∥d w F (w)∥ 2 ). Now, taking in (2.3) the test function ψ among the polynomial basis functions leads to the following linear system allowing the calculation of the solution moments

u c m N k m=1 d u c m dt ωc σ c m σ c p dV = ωc F (u c h ) ∇ x σ c p dV - ∂ωc σ c p F n dS, ∀ p ∈ [1, N k ]. (2.6)
The terms ωc F (u c h ) ∇ x σ c p dV and ∂ωc σ c p F n dS are respectively referred to as volume and surface integrals. In (2.6), we identify ωc σ c m σ c p dV = (M c ) mp as the generic coefficient of the symmetric mass matrix M c . The scheme (2.6) can then be reformulated in a compact matrix-vector form as follows

M c d U c dt = Φ c , (2.7) with (U c ) m = u c
m the solution vector filled with the polynomial moments, and where the so-called DG residuals Φ c write

(Φ c ) m = ωc F (u c h ) ∇ x σ c m dV - ∂ωc σ c m F n dS. (2.8)
Similarly to what has been done in the one-dimensional case in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], let us now demonstrate the equivalency between discontinuous Galerkin schemes and a finite volume like method on a subgrid, and exhibit the corresponding subcell numerical fluxes that will be referred to as high-order reconstructed fluxes. To do so, we first need to subdivide the mesh cells into subcells. Let us emphasize that to obtain a relation of equivalency, one would need the same number of Degrees of Freedom (DoF) as number of subcells. Indeed, it will then be possible to get a one-to-one relation between subcell mean values obtained through a FV-like method provided with reconstructed fluxes and the solution polynomial moments deriving from the DG schemes previously introduced. But the procedure presented in the following is not limited to this case, and this is absolutely possible to exhibit the corresponding reconstructed fluxes even in the case where the number of subcells does not fit the dimension of the polynomial space P k . In the case where more subcells than DoF are used, one can either apply a least square procedure to the set of submean values obtained through the FV-like method to return to a (k + 1) th order DG scheme, or directly return to an higher-order DG scheme in the particular case where the number of subcells corresponds to the dimension of higher degree polynomial space. If one uses less subcells than DoF, then with the same techniques than in the previous case, one can return to a lower DG scheme from the subcell mean values set. In the remainder, for sake of positivity-preserving purpose, we will consider N k subcells.

Even if the choice of the cells subdivision may have an effect of the DG correction technique, for the following theoretical part it has no influence what so ever. The only constrain is that the projection matrix P c further introduced in (2.12) has to been non-singular. Many subdivisions can be found in the literature for other methods relying of a subgrid, as spectral volume methods for instance [START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids: II. Extension to two-dimensional scalar equation[END_REF][START_REF] Harris | Partition design and optimization for high-order spectral volume schemes[END_REF] or subcell shock capturing technique as [START_REF] Huerta | A simple shock-capturing technique for high-order discontinuous Galerkin methods[END_REF]. In Figure 1, three types of subdivision are displayed for a triangular cell, in the case of a 4 th -order DG scheme. Let us emphasize that the one depicted in Figure 1(c) is the most widely used in subgrid techniques, and has the advantages to be invariant by rotation and can be generated for any order of accuracy. For the numerical applications, this subdivision will be compared to a way simpler one, Figure 1(a), which has the benefits to be extremely simple to implement and where the subcells normals are nothing but the primal triangular cell ones.

Even if only triangular grids are considered for numerical applications, let us emphasize that the following demonstration as well as the subcell correction technique presented in this paper are not limited to this case. Any grid made of generic polygonal cells can be considered. Examples of possible subdivisions for DG schemes up to P 3 are displayed in Figure 2. Curvilinear meshes with curvilinear cell subdivision could also be used, and will be the topic of a near future paper. Let us now reformulated DG as a FV-like scheme provided the definition of the so-called reconstructed fluxes.

d u c m dt = 1 |S c m | N k q=1 d u c q dt S c m σ c q dV, (2.10)
which can be put into into a matrix-vector form as follows

d U c dt = P c d U c dt . (2.11)
In (2.11), the projection matrix P c , which has to be invertible for the subdivision to be admissible, is defined

(P c ) mp = 1 |S c m | S c m σ c p dV. (2.12)
Remark 1. Let us note that if the number of subcells does not coincide with the polynomial space dimension N k , then the projection matrix P c is not square. For the subdivision under consideration to be admissible, this projection matrix has to present a left inverse, the reconstruction operator R c , such that R c P c = I N k . In the case of N k subcells, R c reduces to P -1 c .

In relation (2.11)

dt = - 1 |S c m | S v p ∈ V c m f c mp F n dS, (2.15)
where f c mp denotes the face between subcells S c m and S v p . Let us emphasize that S v p ∈ V c m can either be inside cell ω c or in one of its neighbors ω v with v ̸ = c. This situation is displayed in Figure 3 where S c m would be colored red, while its faces neighboring subcells would be colored green. Now, similarly to what has been done in the 1D case, we impose that on the boundary of cell ω c the reconstructed flux coincide with the DG numerical flux

F n| ∂ωc = F n , (2.16)
where F n has been previously introduced in (2.4). This particular choice, while not unique, permits to recover in the most natural way the time evolution governing equation of the cell average value u c

d u c dt := 1 |ω c | ωc ∂ u c h ∂t dV = - 1 |ω c | ∂ωc F n dS,
obtained through DG scheme definition (2.3) with ψ = 1. By means of (2.16), definition (2.15) then rewrites

d u c m dt = - 1 |S c m |    S c p ∈ ( V c m f c mp F n dS + ∂S c m ∩∂ωc F n dS    , (2.17)
where ( V c m stands for the set containing only the face neighboring subcells of S c m inside ω c . For now on, an orientation will be assigned to each face. Then, taking two subcells S c m and S v p , we introduce the sign function ε c mp depending on the orientation of face f c mp

ε c mp =          1 if face f c mp is direct or if f c mp ⊂ ∂ω c , -1 if face f c mp is indirect, 0 if S v p / ∈ V c m .
(2.18)

Obviously, ∀S c p ∈ ( V c m , we have ε c pm = -ε c mp . Now, unlike the 1D case where a pointwise definition of the reconstructed flux was given inside the cell, we make use here of a face integrated version of the high-order DG reconstructed flux. Indeed, for a face f c mp , let F mp be defined as follows

f c mp F n dS = ε c mp F mp . (2.19)
Since the face orientation has been carried through ε c mp , the face integrated quantity F mp is then continuous, i.e. ∀S c p ∈ ( V c m , F pm = F mp . Now, denoting by N c f the number of subcells faces inside ω c , meaning not belonging to ∂ω c , let us introduce F c ∈ R N c f the vector containing all the interior faces reconstructed fluxes. The subcell mean values governing equations (2.17) then yield the following system

-A c F c = D c d U c dt + B c , (2.
-A c F c = Q c M -1 c Φ c + B c . (2.21)
To solve such system, and obtain an explicit expression of the reconstructed fluxes F c , we make use of the technique employed in [START_REF] Abgrall | Some Remarks About Conservation for Residual Distribution Schemes[END_REF] in a similar context. Since the graph is connected, the rank and the kernel of the Laplacian matrix of the interior subgrid graph L c = A c A t c respectively are N k -1 and span{1}. Then the inverse of L c on the orthogonal of its kernel is defined, and is henceforth denoted by L -1 c . For any λ ̸ = 0, this inverse writes

L -1 c = (L c + λ Π) -1 - 1 λ Π, (2.22) with Π = 1 N k (1 ⊗ 1) ∈ M N k .
Finally, we are able to exhibit the following definition of the reconstructed flux

F c = -A t c L -1 c Q c M -1 c Φ c + B c . (2.23)
Remark 2. In the definition of the reconstructed flux through DG residual, (2.23), the only terms depending on time are Φ c , the residual which is directly available in any DG code, and the boundary contribution B c . All the other terms can be evaluated initially, once and for all. Furthermore, let us emphasize that if all the mesh cells have the same structure, as triangles for example, then by means of a mapping to a referential triangle, one can prove that the projection matrix P c , the adjacency matrix A c and the generalized inverse of the graph Laplacian matrix L -1 c do not depend on the cell under consideration, but only on the order of approximation and on the choice of the subdivision. Remark 3. As said previously, for the subdivision to be admissible, the projection matrix has to be invertible, which is almost always true. However, the type of cell subdivision may have an influence on the condition number of the projection matrix, and then may introduce more numerical errors during the simulation. In Section 4, the projection matrix condition numbers are gathered for the four types of cell subdivisions used in practical applications, for different orders of accuracy, see Table 2.

Once we have computed the reconstructed flux F c , we can simply recover the polynomial solution governing equation as follows

d U c dt = -Q -1 c A c F c + B c . (2.24)
Let us gather all the previous results into the following theorem.

Theorem 2.1. Discontinuous Galerkin schemes, expressed in the following equation in cell ω c in a compact matrix-vector form through residuals

d U c dt = M -1 c Φ c , (2.25)
where the DG residual Φ c has been defined in (2.8), can be recast into N k subcell finite volume like schemes as

d U c dt = -D -1 c A c F c + B c , (2.26)
where the finite volume fluxes F c , referred to as reconstructed fluxes, is defined through equation (2.23). All the other quantities involved in (2.25) and (2.26) have been previously introduced. Now, for a deeper understanding of the reconstructed flux F c , let us seek its relation with the interior flux F (u c h ) and the DG numerical flux F n . 2.2. Reconstructed flux through fluxes. In this section, we aim at expressing the reconstructing fluxes only through the interior flux F (u c h ) and a correction term taking into account the jump at the cell boundaries, similarly to what is done in SBP operator with SAT boundary treatment [START_REF] Fisher | High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains[END_REF][START_REF] Gassner | A Skew-Symmetric Discontinuous Galerkin Spectral Element Discretization and Its Relation to SBP-SAT Finite Difference Methods[END_REF] or CPR schemes [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF][START_REF] Wang | A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids[END_REF]. To do so, we will not make use of the definition of DG scheme through residual (2.25) but through fluxes (2.3). The first step is to substitute in (2.3) the interior flux F (u c h ) by F c h its L 2 projection onto P k+1 (ω c )

2 . If one uses nodal DG or any collocation of the interior flux, this step can obviously be skipped.

Performing a second integration by parts leads to the so-called strong form of DG scheme

ωc ∂ u c h ∂t ψ dV = - ωc ψ ∇ x F c h dV + ∂ωc ψ (F c h n -F n ) dS, ∀ ψ ∈ P k (ω c ). (2.27)
Similarly to what has been done in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] for the one-dimensional case, let us introduce the N k sub-resolution basis functions {ϕ m } m . These particular basis functions of P k (ω c ), which can be seen as the L 2 projection of the subcell indicator functions

1 S c m (x) onto P k (ω c ), are defined such that ∀ ψ ∈ P k (ω c ) ωc ϕ m ψ dV = S c m ψ dV, ∀ m = 1, . . . , N k . (2.28)
Now, because equation (2.27) holds for any polynomial function ψ of degree k, let us substitute ϕ m for ψ in DG schemes

ωc ∂ u c h ∂t ϕ m dV = - ωc ϕ m ∇ x F c h dV + ∂ωc ϕ m (F c h n -F n ) dS. (2.29)
Through the sub-resolution property (2.28), equation (2.29) can be recast as

|S c m | d u c m dt = - S c m ∇ x F c h dV + ∂ωc ϕ m (F c h n -F n ) dS. (2.30)
Finally, we obtain the following governing equation of the submean values time evolution

d u c m dt = - 1 |S c m | ∂S c m F c h n dS - ∂ωc ϕ m (F c h n -F n ) dS . (2.31)
The use of the reconstructed flux definition (2.14) directly leads to the following relation between the different types of fluxes

∂S c m F n dS = ∂S c m F c h n dS - ∂ωc ϕ m (F c h n -F n ) dS. (2.32)
One can see how the reconstructed flux is connected to the interior polynomial flux and the jump at the cell interface between the interior flux and the DG numerical flux. To further develop relation (2.32), we make the same assumption on the cell boundary reconstructed fluxes as before (2.16). Expression (2.32) then rewrites

∂S c m \∂ωc F n dS = ∂S c m \∂ωc F c h n dS - ∂ωc ϕ m (F c h n -F n ) dS, (2.33)
where, for sake of conciseness, ϕ m reads as follows

ϕ m = ϕ m if x ∈ ∂ω c \ ∂S c m , ϕ m -1 if x ∈ ∂ω c ∩ ∂S c m .
(2.34)

Similarly to (2.19), let F mp be the face integrated integrated value of the polynomial interior flux

f c mp F c h n dS = ε c mp F mp . (2.35)
Then, if F c is the vector containing all the interior faces fluxes, one gets

A c F c = A c F c -G c , (2.36)
where the G c contains the boundary contribution as

(G c ) m = ∂ωc ϕ m (F c h n -F n ) dS. (2.37)
Finally, by means of the same graph Laplacian technique used previously, we are able to express the reconstructed flux through the interior flux and a boundary correction term

F c = F c -A t c L -1 c G c . (2.38) Remark 4. We can rewrite (2.38) as F c = F c -E(F c h n -F n )
, where E(.) would be a correction function taking into account the jump between the polynomial flux and the numerical flux on the cell boundary. This permits to demonstrate once more that in DG schemes the numerical diffusion deriving from the jump in term of flux across the cell interfaces is distributing elsewhere inside the cell, here at the subcells faces. The sub-resolution basis functions act as weighted functions in the diffusion distribution. Let us note that another choice in the correction term function E(.) would lead to a different scheme. For instance, setting E(.) = 0 leads to spectral volume scheme of Z.J. Wang [START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids: Basic Formulation[END_REF].

As previously, let us gather the previous results in another theorem.

Theorem 2.2. Discontinuous Galerkin schemes, expressed in the following equation in cell ω c through volume and boundary flux contribution

ωc ∂ u c h ∂t ψ dV = ωc F (u c h ) ∇ x ψ dV - ∂ωc ψ F n dS, ∀ ψ ∈ P k (ω c ), (2.39)
can be recast into N k subcell finite volume like schemes as

d u c m dt = - 1 |S c m | S v p ∈ V c m ε c mp F mp , (2.40)
where the finite volume fluxes F mp , referred to as reconstructed fluxes, are either defined in equation (2.38) if S v p ⊂ ω c and in (2.16) otherwise. Both theorems 2.1 and 2.2 are perfectly equivalent. While the definition of DG reconstructed flux (2.38) enables a better comprehension of how those fluxes are related to the interior flux F (u c h ) and how the numerical diffusion deriving from the jump at the interface is distributed inside the cell, the definition of reconstructed fluxes through residual (2.23) is a lot easier to implement and do not require the definition of sub-resolution basis functions.

A posteriori local subcell correction.

The previous reformulations of DG scheme into subcell FV-like scheme through the definition of reconstructed fluxes enable us to construct our a posteriori local subcell correction (APLSC). In few words, the reconstructed fluxes F mp will be modified in a robust way in subcells where the original DG scheme has failed. Let us emphasize that the popularity and the number of subcell correction techniques have extensively grown these past years, see [START_REF] Huerta | A simple shock-capturing technique for high-order discontinuous Galerkin methods[END_REF][START_REF] Sonntag | Shock capturing for discontinuous Galerkin methods using finite volume subcells[END_REF][START_REF] Dumbser | A simple robust and accurate a posteriorisub-cell finite volume limiter for the discontinu-ousGalerkin method on unstructured meshes[END_REF][START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF][START_REF] Kuzmin | Subcell flux limiting for high-order bernstein finite element discretizations of scalar hyperbolic conservation laws[END_REF][START_REF] Pazner | Sparse invariant domain preserving discontinuous galerkin methods with subcell convex limiting[END_REF][START_REF] Hennemann | A provably entropy stable subcell shock capturing approach for high order split form dg for the compressible euler equations[END_REF]. Those shock capturing and property preserving methods generally rely on a low-order scheme combined, at the subcell level, with a high-order one. However, let us note that in most cases, all the subcells contained in a cell will be impacted if something bad happened somewhere in the cell. In [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], we have introduce, for the one-dimensional case, a new technique permitting to correct the solution in a subcell without modifying the solution elsewhere. This particular feature allow to retain as much as possible the very precise subcell resolution of high-order DG schemes. The present paper aims at presenting the two-dimensional version of this correction. Let us emphasize that, up to our knowledge, this is the only technique working on totally unstructured grids and permitting the modification of the scheme, locally at the subcell level, without impacting the solution everywhere in the cell under consideration.

Let us mention that until now, only the semi-discrete version of schemes and their corresponding analysis were presented. To achieve high-accuracy in time, we make use of SSP Runge-Kutta time integration method [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. But, in the light of the fact that these multistage time integration methods write as convex combinations of first-order forward Euler scheme, the correction DG procedure will be presented for the simple case of this latter time numerical scheme, for sake of simplicity. DG schemes (2.3) provided with first-order forward Euler time integration writes,

∀ ψ ∈ P k (ω c ) ωc u c, n+1 h ψ dV = ωc u c, n h ψ dV + ∆t ωc F (u c, n h ) ∇ x ψ dV - ∂ωc ψ F n dS . (3.1)
The numerical solution u c, n h on cell ω c , being assumed to be a P k polynomial, is uniquely defined through its mean values on the N k subcells, and reversely. To go from one representation to another, namely from {u c, n m } m the solution moments to {u c, n m } m the subcell mean values, one has just to make use of the projection matrix P c or its inverse

P c    u c, n 1 . . . u c, n N k    =    u c, n 1 . . . u c, n N k    . (3.2)
Let us now introduce the correction procedure. First, we assume that at time t n the numerical solution u n h is satisfactory in the sens that, on any cell ω c , the subcell mean values are admissible regarding some criteria yet to be defined. Then, we compute u n+1 h a candidate solution through the uncorrected DG scheme. The third step is then crucial. Indeed, we then have to check if the new uncorrected solution is admissible. If it is the case, we can go further in time without any special treatment. Otherwise we have to return to time t n and recompute the solution locally by means of a more robust scheme. This step is crucial in the sens that it will tell us if and where a new computation would be required.

Troubled zone detector.

Regarding the troubled zone detectors, we simply extend to the twodimensional unstructured case the ones used in the one-dimensional frame [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. In this previous work, two detection criteria were mainly used, namely one ensuring the physical admissibility of the numerical solution (PAD) and another addressing the apparition of spurious oscillations (NAD). Let us then recall these two criteria.

Physical admissibility detection (PAD).

• Check if the different mean values u c, n+1 m lie in a chosen convex physical admissible set (maximum principle for SCL, positivity of the pressure and density for Euler, . . . ). Entropy stability can even be added to this admissible set.

• Check if there is any N aN values Those are the minimum requirements if one wants to enforce code robustness. Now, in order to tackle the issue of spurious oscillations, we make use of a local maximum principle. Indeed, through the respect of the CFL, the solution in cell ω c at time t n+1 has to remain in the bounds of the solution at the previous time step t n wherein ω c and its first face neighbors. This condition is reformulated in the following detection criterion.

Numerical admissibility detection (NAD).

• Check if the following Discrete Maximum Principle (DMP) on submean values is ensured:

min v ∈ N (S c m ) (u n v ) -δ ≤ u c,n+1 m ≤ max v ∈ N (S c m ) (u n v ) + δ,
where N (S c m ) is some set of S c m neighboring subcells, including subcell S c m . The parameter δ, although not required, enables a small relaxation of the NAD criterion. Similarly to [START_REF] Dumbser | A simple robust and accurate a posteriorisub-cell finite volume limiter for the discontinu-ousGalerkin method on unstructured meshes[END_REF], we set δ as follows

δ = max 10 -4 , 10 -3 max v ∈ N (S c m ) (u n v ) -min v ∈ N (S c m ) (u n v ) . (3.3)
Remark 5. The bigger the set N (S c m ) is, the softer the NAD criterion will be. Reversely, a smaller set would lead to more subcells considered as problematic. Different sets will be considered for linear and non-linear problems, see Sections 4.1 and 4.2.

Let us enlighten that NAD criterion relies on a maximum principle based on subcell mean values. And because this maximum principle is not concerned with the whole polynomial set of values, it is very wellknown that one has to relax it to preserve scheme accuracy in the presence of smooth extrema. To do so, we make use of the two-dimensional version of the one introduced in the 1D case, [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF].

Detection of smooth extrema. This smooth extrema detection criterion is based on an idea present in different limitations, as the hierarchical slope limiter [START_REF] Kuzmin | A vertex-based hierarchical slope limiter for p-adaptative discontinuous Galerkin methods[END_REF]. In this limitation technique, the numerical solution is supposed to exhibit a smooth extrema if at least the linearized version of the numerical solution spatial derivatives, i.e.

(3.4a) (3.4b)      v c x (x) = ∂ x u c, n+1 h + ∇ x (∂ x u c, n+1 h ) (x -x c ), v c y (x) = ∂ y u c, n+1 h + ∇ x (∂ y u c, n+1 h ) (x -x c ),
present a monotonous profile. In (3.4), x c denotes the centroid of cell ω c , while ∂ x\y u c, n+1 h and ∇ x (∂ x\y u c, n+1 h ) are nothing but the averaged values on ω c of the successive partial derivatives of u c h . In practice, the NAD relaxation used here works as a vertex-based limiter on v c x\y (x). Due to their linearity, functions v c x\y (x) attain their extrema at the vertices x p ∈ P c , where P c stands for the set of vertices of cell ω c . Then, we consider that the exact weak solution underlying the numerical solution u h presents a smooth profile in cell ω c if, for any vertex x p ∈ P c , the linearized spatial derivative functions ensure the following constraint In (3.6), V p represents the set of cells that share x p as a vertex, i.e. ω v ∈ V p =⇒ x p ∈ P v . Practically, if for any vertex x p ∈ P c , conditions (3.5) are ensured, we then consider that the numerical solution presents a smooth profile on cell ω c . In this particular case, the NAD criterion is relaxed allowing the preservation of smooth extrema along with the order of accuracy for smooth problems, see Section 4. Remark 6. We have presented here the detection based on the linearized first-order derivatives of the solution. This would work for any higher order derivative. Furthermore, such relaxation procedure can also be applied at the subcell level to also preserve smooth extrema even within a cell at a smaller length scale. This would be useful in the context of coarse grids. Actually, because the subcell smooth extrema relaxation technique works well for both coarse and fine meshes, this will be the procedure used for the numerical applications, Section 4. Now that we have detailed the troubled subcell detector, the correction procedure will be presented in the next subsection.

v min x, p ≤ v c x (x p ) ≤ v max x, p and v min y, p ≤ v c y (x p ) ≤ v max y, p , (3.5
3.2. Correction. In this subsection, the original correction procedure introduced in the one-dimensional frame in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] will firstly be extended to the two-dimensional unstructured grid case. Then, another technique blending in a convex combination fashion DG reconstructed fluxes and first-order FV fluxes for admissible subcells in the vicinity of troubled areas will also be presented and prove to improve results in comparison to the original algorithm.

Original correction principle.

The very simple idea that forms the basis of the original correction procedure introduced in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] is the following: if the uncorrected DG scheme has produced a numerical solution u c, n+1 h on cell ω c , which is not admissible in subcell S c m in regards to the detection criteria presented previously, the subcell mean value u c, n+1 m will be recomputed by means of a more robust scheme. To do so, and because uncorrected DG scheme is equivalent to subcell finite volume scheme with the appropriate high-order reconstructed fluxes, see Theorems 2.1 and 2.2, we substitute on the boundaries of subcell S c m the high-order reconstructed fluxes with some first-order finite volume numerical fluxes. The submean value u c, n+1 m will then be recomputed by means of a simple and robust first-order finite volume scheme. This concept is depicted in Figure 3, where the troubled subcells are colored red. Because it is of fundamental importance to preserve scheme conservation, the first face neighboring subcells, colored green in Figure 3, have to be also recomputed since we have modify the reconstructed fluxes on the boundary of the troubled subcell. The submean values of the neighboring subcells are then computed through a FV-like scheme with first-order numerical flux on one or more faces and high-order reconstructed fluxes on the remaining interfaces. For the remaining subcells, colored gray in Figure 3, because the corresponding reconstructed fluxes have not been modified, there is no need to recompute them. The corresponding submean values are hence the values obtained through the uncorrected DG scheme. This original concept of the a posteriori local subcell correction of DG scheme is summarized in the following flowchart:

[lef t = {] F mp = ε c mp l c mp F u c,n m , u v,n p , n mp if S c m or S v p ∈ V c m is either marked, F mp = F mp otherwise,
where l c mp is the length of face f c mp .

5. Through the corrected reconstructed flux, recompute the submean values for tagged subcells and their first neighboring subcells as In light of this correction procedure flowchart, it is clear that the DG solution will only be affected at the subcell scale. Furthermore, the corrected scheme is conservative at the subcell level by construction.

u c, n+1 m = u c, n m - ∆t |S c m | S v p ∈ V c m ε c mp F mp .
Remark 7. In the present correction, when needed we substitute the reconstructed flux value F mp with a first-order numerical flux F u c,n m , u v,n p , n mp . Obviously, other choices are possible and may even be more appropriate, as for instance a second-order TVD numerical flux or even a WENO numerical flux. In the 1D, results with both first-order and second-order TVD corrections were presented, see [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. In the present article, for sake of simplicity and conciseness, only the first-order correction previously detailed will be used.

Remark 8. Let us emphasize that since this a posteriori correction is based on first-order finite volume scheme, maximum principle or positivity preservation in the case of systems are enforced by construction. However, it is crucial to note that to start from a discrete representation of the initial datum that respects its bounds (or ensures positivity), the initialization has to be carried out through the integration of the initial solution on the subcells to obtain the different submean values, and then through relation (3.2) compute the corresponding polynomial moments of the solution. Traditionally, in DG schemes the initialization is handled by either assigning at the solution points the initial datum value to the numerical solution, or by a L 2 projection onto P k . In either of these procedures, nothing ensures that the submean values would respect the bounds of the initial datum.

Let us recall that in the 1D case [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], for non-linear problems using very high-order schemes on coarse grids, the numerical solution proved to remain slightly oscillatory at the subcell level. To overcome this issue, we were artificially enlarging the stencil to correct by also marking, for k ≥ 3, the first face neighboring subcells of a troubled subcell. Now, in this 2D frame, a modified algorithm making use of convex combination of DG reconstructed fluxes and first-order FV fluxes for admissible subcells in the vicinity of troubled areas will be presented and prove to produce better results than the original correction, in this context of very high-orders and coarse meshes.

New correction principle.

As said previously, for very high-order DG methods on coarse grids, the corrected scheme may locally jump, in two subcells scale, from a very precise approximation to a robust but very low accurate first-order representation. In the particular situation, numerical experiments have proved that some subcell oscillations may remain due to this phenomenon. To overcome this issue, we will now introduce a new correction algorithm in which the stiffness of the order transition will be smoothened by means of convex combination between DG reconstructed fluxes and first-order FV fluxes as follows

F mp = θ mp ε c mp l c mp F u c,n m , u v,n p , n mp + (1 -θ mp ) F mp , (3.7)
where θ mp is a function of the distance to a non-admissible subcell. Obviously, the farther from the troubled subcell, we are the less of first-order FV we want. Many smoothness indicators exist in the literature, see for instance [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V-A second-order sequel to Godunov's method[END_REF][START_REF] Harten | ENO schemes with subcell resolution[END_REF][START_REF] Burbeau | A problem-independent limiter for high-order Runge Kutta discontinuous Galerkin methods[END_REF][START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF][START_REF] Persson | Sub-cell shock capturing for discontinuous Galerkin methods[END_REF][START_REF] Balsara | A sub-cell based indicator for troubled zones in RKDG schemes and a novel class of hybrid RKDG+HWENO schemes[END_REF], and could be used in this context to determine a relevant bending coefficient. Another way around could be to adopt flux limiting or Flux-Corrected transport (FCT) approach, see [START_REF] Boris | Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works[END_REF][START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF][START_REF] Kuzmin | Flux-Corrected Transport[END_REF][START_REF] Kuzmin | Flux-Corrected Transport: Principles, Anlgorithms and Applications[END_REF][START_REF] Lohmann | Flux-corrected transport algorithms for continuous Galerkin methods based on high order Bernstein finite elements[END_REF][START_REF] Kuzmin | A new perspective on flux and slope limiting in discontinuous Galerkin methods for hyperbolic conservation laws[END_REF][START_REF] Kuzmin | Bound-preserving flux limiting for high-order explicit Runge Kutta time discretizations of hyperbolic conservation laws[END_REF], to find the proper blending ensuring a discrete maximum principle. To remain as simple as possible, and because this paper is concerned with the introduction of the DG-FV equivalency and the basic correction principle on 2D unstructured grids, we make use here of a very naive procedure. The principle is the following, for marked subcells detected through the troubled subcell indicator, first-order FV numerical flux is used on its boundaries, i.e. θ mp = 1 in (3.7). Then, for its face neighboring subcells S v p ∈ V c m , convex combination (3.7) with θ mp = 3 4 would be chosen for their remaining boundaries. Now, introducing V c m the set containing both the face and node neighboring subcells, fluxes on faces of S v p ∈ V c m \ V c m which yet have not been corrected will be defined through at the blending coefficient θ mp = 1 2 . Finally, the remaining face neighbors of S v p ∈ V c m will see their remaining boundaries associated with a numerical flux calculated through (3.7) with θ mp = 1 4 . This naive technique is displayed in Figure 4. One could think that doing so, we have substantially enlarged the stencil of subcells to be corrected, and thus reduce the simulation code efficiency. However, in practice it is generally not the case, as one can see in Figure 6. Moreover, the computational time will also be generally slightly reduced. Emphasizing that the correction is done in an a posteriori fashion and thus has to be potentially iterated multiple times at a time step to reach an admissible solution (see point 7. of the previous flowchart), if the stiff original correction introduces small oscillations at the subcell scale, the stencil will be automatically enlarged along with an increase of the computational time, through the correction iteration. With the new correction principle, only one iteration is generally needed. So even if at first, the set containing marked subcells is smaller in the original correction, it will end up with approximately the same size than with this new approach, and required more iterations and thus more computational effort. To assess this matter and compare both corrections, we make use of the Burgers equation, defined through (2.1a) and flux function F (u) = 1 2 u 2 , u 2 t , with the smooth initial solution u 0 (x) = sin(2π (x + y)). The domain is chosen as the unit square [0, 1] 2 with periodic boundary condition. Through time, the exact solution will exhibit two stationary shocks along the lines defined by x ∈ [0, 1] 2 , x + y = 0.5 and x ∈ [0, 1] 2 , x + y = 1.5 , see Figure 5.

F c F c 3 4 F c + 1 4 F c 1 2 F c + 1 2 F c 1 4 F c + 3 4 F c (
We run this test case with a sixth order DG scheme, on a quite coarse unstructured triangular grid made of 576 cells, with both a posteriori local subcell corrections. We make use here of the simple structured subdivision introduced in Figure 1(a). To compare the two approaches, let us first display the marked subcells to be corrected. In Figures 6 and7, we color all the subcells that have been corrected in the different Runge-Kutta steps during the last time step. Firstly, one can note that both corrections, through the NAD criterion, have accurately capture the discontinuities, as the subcells to be corrected remain in a small vicinity of the shocks. Secondly, we can also observe how these corrections operate locally at the subcell scale. Finally, as depicted in Figure 6, the number of subcells to be corrected remains approximately the same with both approaches. Actually, with the original correction on average 10% of the total number of subcells have to be corrected through the computation, while 14% of the subcells with the new approach. But, even if indeed the number of subcells to be corrected has slightly increases through this new procedure, if we compare computational efficiency, it took 2 minutes and 22 seconds to the code with the original correction, and 2 minutes 15 seconds with the new one. As said previously, it comes from the fact that more correction iterations are required with the original procedure. For this calculation, on average 2.86 iterations are needed with the original approach when the correction has been triggered, for a maximum of 6 iterations, while only 1.46 iterations with the new approach, with a maximum of 3 iterations during the whole calculation. For a better understanding of how work those two different approaches, we zoom on the area [0.65, 0.9] 2 . Let us first note that because the dimension of P 5 is 20, each cell has been accordingly subdivided into 20 subcells. In the light of Figure 7, one can observe that a lot less subcells are corrected through a purely first-order finite volume with the new correction than with the original one, which enables even more the preservation of the high accurate subcell resolution of high-order DG schemes. In Figure 7(a), one can see the two types of subcells to be corrected, meaning the troubled subcells colored red and their face neighbors colored green which have to be recalculated to preserve the scheme conservation. In Figure 7(b), there are five types of subcells to be corrected, from the troubled subcells to be recomputed through a first-order FV scheme to the magenta, purple and blue ones which are recalculated through a convex combination between first-order FV numerical flux and the high-order DG reconstructed flux, with different weights. The green ones are the face neighbors of the previous subcells to be recomputed to preserve scheme conservation. While the previous results showed how those two approaches work and do have a quasi-equivalent computation cost, let us now compare the approximated solutions obtained and assess the benefit of the new correction. In Figure 8, the subcell mean values obtained by means of a sixth-order DG scheme corrected through the original and the new procedures are displayed. Comparing Figures 8(a) and 8(b), we can see that the new correction has improved the accuracy of the scheme by a sharper and less oscillatory representation of the shock, despite the very coarseness of the mesh used.

Finally, we plot in Figure 9 for both corrections, as well as for the exact entropic solution, the subcell mean values u c m versus r c m = x c + y c -1, where (x c , y c ) stands for the barycenter of subcell S c m . These results demonstrates once more how the new correction leads to a better representation of the shocks, with a less oscillatory profile. For that reason, the new correction based on a convex combination of first-order FV numerical flux and high-order DG reconstructed flux, with decreasing weight in the vicinity of a troubled subcell, will be adopted for the numerical applications presented in the next section. 4. Numerical results. In this numerical results section, we make use of several widely addressed and challenging test cases to demonstrate the performance and robustness of this a posteriori local subcell correction for discontinuous Galerkin schemes. In all following test cases, the simple case of global Lax-Friedrichs numerical flux will be used for both the DG scheme and the first-order finite volume reconstructed flux correction. Regarding the cell decomposition into subcells, while this has no impact on the reformulation of DG schemes into subcell finite volume method, see Section 2, it may have a slight impact on the results obtained by means of the present correction. In the one-dimensional case, see [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], it has been demonstrated that the use of a non-uniform subdivision, for instance by means of the Gauss-Lobatto points, leads to better results compared to a uniform subdivision. In this two-dimensional framework the two types of subdivisions introduced respectively in Figure 1(a) and 1(c).

Regarding the time integration, we make use of the classical third-order SSP Runge-Kutta scheme, see for instance [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. As the correction described earlier combines both DG scheme on the primal cells ω c and FV scheme on the subcells S c m , the time step is computed adaptively using the following CFL condition We recall that l c mp stands for the length of face f c mp separating subcell S c m and its neighbor S v p , while l cv is the length of the interface between cell ω c and its neighbor ω v . Let us note that in cases where we compute rates of convergence, a time step ∆t ≤ min

c d k+1 3 c
is used in order to make the time error negligible in comparison to the spatial discretization error.

Let us emphasize that in all figures to come, if not stated otherwise, the solution subcell mean values are displayed. Consequently, for a mesh made of N c cells, the numerical solution will be represented on N k N c subcells, where N k = (k+1)(k+2) 2 .

Linear case. Let us first assess the performance and accuracy of the present A Posteriori Local

Subcell Corrected DG (APLSC-DG) scheme in the case of 2D linear conservation laws. In Section 3.1, when the NAD criterion based on a discrete maximum principle has been introduced, we did not specified the set N (S c m ) of subcell S c m characterizing the DMP. For linear problems, we make use of a cell-wise DMP, meaning N (S c m ) will be constituted by all the subcells of cell ω c , as well as the subcells of ω v the face neighboring cells of ω c . By means of notations previously introduced, this definition can be rewritten as

N (S c m ) = S v q ; ∀ ω v ∈ V c ∪ {ω c }, ∀ q ∈ 1, N k . (4.3)
This particular set is depicted in Figure 10, for both the simple structured subdivision as well as the polygonal Voronoi-type one. In this figures, the subcell S c m under consideration would be colored red, while the subcells constituting N (S c m ) would be colored green. Let us emphasize that subcell S c m is also part of N (S c m ). 4.1.1. Linear advection. To display the efficiency of DG schemes plus correction, let us first assess how the APLSC behaves in the linear advection case. To this end, we consider the following equation

(4.4a) (4.4b)    ∂ t u(x, t) + A(x) ∇ x u(x, t) = 0, (x, t) ∈ [0, 1] 2 × [0, T ], u(x, 0) = u 0 (x), x ∈ [0, 1] 2 ,
where the transport velocity will be set to A = (1, 1) t .

Linear advection of a smooth signal. We start from a smooth initial datum u 0 (x, y) = sin(2π (x+y)), and consider periodic boundary conditions. We assess the scheme accuracy after one period, namely at time t = 1. In Figure 11, the numerical solution of the APLSC of sixth-order DG scheme, obtained on grid made of only 100 cells, is depicted. One can see that with only 100 cells, the corrected DG scheme is extremely accurate. Actually, the correction procedure does not activate in this case, which proves that the relaxation criterion on smooth extrema works properly. The rates of convergence are gathered in Table 1 and do exhibit a convergence to six.

L 1 L 2 L ∞ h E h L1 q h L1 E h L2 q h L2 E h L∞ q h L∞ 1 10
1.62E-7 6.00 1.81E-7 6.00 3.98E-7 5.96

1 20
2.53E-9 5.97 2.82E-9 5.96 6.38E-9 5.41 1 Convergence rates for the linear advection case for a 6th-order APLSC-DG scheme Linear advection of a crenel signal. To assess the efficiency of the correction presented in the presence of discontinuity, let us start with the simple case of the advection of a crenel signal, where the initial solution u 0 is defined as follows

u 0 (x) =          1 if (x + y) ∈ [ 1 4 , 1 2 ] ∪ [ 5 4 , 3 2 ], 0 if (x + y) ∈ [ 3 4 , 1] ∪ [ 7 4 , 2], 1 2
otherwise.

(4.5)

In Figure 12, we compare the uncorrected and corrected versions of the 6th-order DG scheme on an unstructured grid made of 576 cells, after one period. One can see that in both cases, the numerical solution obtained is very accurate, but the one obtained through the APLSC-DG method respects the maximum principle as the final solution remains in the bounds of the initial one. In Figure 13, the submean values versus (x + y -1) coordinate are compared for both the uncorrected and the APLSC-DG schemes. We note how the spurious oscillations in the vicinity of the discontinuities have been removed, while preserving the very precise resolution of the 6th-order DG scheme.

In the previous examples, the simple uniform structured cell subdivision has been used. Let us emphasize that the cell subdivision does not theoretically impact the equivalency between DG scheme and subcell FVlike scheme defined through the reconstructed fluxes introduced previously, see Theorems 2.1 and 2.2. This choice, if it does, should only influence the corrected scheme. To illustrate such statement, let us run the linear advection of the crenel signal on five periods for different types of cell subdivision.

In Figure 14, the simple case of structured subdivision is depicted where the cell is partitioned as a quadrilateral cell would be and then split into two to fit the triangular shape. In this case, one can note that the subdivision is not rotation invariant and that a choice has to be made as one corner subcell is a quadrilateral while the other two are triangles. In this work, we made the choice to start the structured cell subdivision from the wider angle corner, which induces the quadrilateral subcell to stand as this particular corner. In Figure 14(a), the cell boundary points defining the subdivision are distributed in an uniform manner, while in Figure 14(b) there are defined as Gauss-Lobatto quadrature points. Let us emphasize that in both cases, those subdivisions are very simple to implement and to generalize to any order of accuracy. Furthermore, the subcells normals are nothing but the ones to the original triangular cell, and the subcells are either triangles or parallelograms. In Figure 15, polygonal subdivisions are displayed. Those Voronoi type subdivisions are widely used in subcell techniques and can be found for instance in [START_REF] Chen | Partitions of a simplex leading to accurate spectral (finite) volume reconstruction[END_REF][START_REF] Huerta | A simple shock-capturing technique for high-order discontinuous Galerkin methods[END_REF] and references within. In the four calculations lead to exactly the same numerical solution.

That being said, the type of subdivision does have an impact on the condition number of the projection matrix, which may lead for very high-order of accuracy to greater numerical errors. In Table 2, the projection matrix condition numbers for the four different subdivisions displayed before and for different orders of approximations are gathered. These different condition numbers have been computed through the ||.|| ∞ consistent norm, such that ||P || ∞ = max i j P ij . It appears that, even if the condition number are closed for these subdivisions, the structured ones do produce higher condition number than the polygonal ones, as well as the uniform ones compared to the non-uniform ones. One can see that higher orders of approximation yield greater condition number. This may ultimately lead to potential issue in the case of extremely high-orders, but in all the cases studied nothing symptomatic has been observed.

P 0 P 1 P 2 P 3 Unif. struct.
While the cell subdivision does not affect the uncorrected DG numerical solution, see Figure 18, it has been proved in the 1D case [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] that it does have an impact on the quality of the results when the correction is used, especially in the linear advection case. To assess if any similar phenomenon exists in the 2D unstructured case, we use the same setup as before but this time with the full APLSC-DG method, see Figures 19,20 and 21. In the light of those results, it appears that uniform subdivisions lead to a way better resolution of the problem under consideration. This huge difference mainly derives from the application of the NAD troubled subcell detector, as if we make just use of the PAD criterion to only enforce the global maximum principle, the difference in the results is a lot more slight. NAD criterion based on a discrete maximum principle is here highly dependent of the subcell aspect ratio. Similar conclusion can be obtained considering other type of problem, even if the difference in the quality of the results would be a lot less significative.

4.1.2. Solid body rotation. We make use of the classical test case taken from [START_REF] Leveque | High-resolution conservative algorithms for advection in compressible flow[END_REF]. Let us then consider (4.4) with a divergence-free velocity field corresponding to a rigid rotation, defined by A(x) = ( 12 -y, x -1 2 ) t . We apply this solid body rotation to the initial data displayed in Figure 22(a), which includes both a plotted disk, a cone and a smooth hump. In Figure 22, the submean values obtained through the 6th-order APLSC-DG scheme on a 576 cells coarse grid are displayed. The uniform structured cell subdivision has been used. One can see on Figure 22(b) how the corrected DG scheme produces a very accurate solution, even using a quite coarse mesh, while still ensuring a global maximum principle as well as a mainly non-oscillatory behavior. In Figure 23, cross-sections of the solution along lines y = 0.25 and x = 0.25 have been plotted. Those results further demonstrates the very high capability of the correction procedure presented. Now, similarly to the linear advection equation, we want to understand how the choice in the cell subdivision impacts the quality of the results now considering a solid body rotation problem. To do so, we compare the solutions obtained by means of the 4th-order APLSC-DG scheme and the four different subdivisions introduced before, see Figures 14 and15, after five full rotations. In the light of Figures 24 and25, we once again note how more uniform cell subdivision lead to less diffused solutions. The solutions cross-sections displayed in Figure 26 confirm this statement. However, it is worth mentioning that the difference in the result is not as critical as it was in the linear advection case. Furthermore, even if the uniform structured subdivision is not rotation invariant, it led to comparable results to the uniform Voronoi-type subdivision, while being a lot more simpler to implement and to generalize to any order of accuracy.

Non-linear case.

Let us now assess the performance and accuracy of our APLSC-DG technique in the 2D non-linear case. Both cases of scalar conservation laws as well as system of conservation laws will be addressed. Similarly to the subsection devoted to the linear case, let us defined set N (S c m ) when considering the NAD criterion on subcell S c m . Unlike the linear case, we make use here of a subcell-wise DMP, meaning N (S c m ) will be constituted by subcell S c m , as well as all its face and node neighboring subcells of S v q , either they belong to the same cell or not. By introducing P c m the set of vertices of subcell S c m as well as N p the set of subcells that share x p as a vertex, this definition can be rewritten as This particular set is depicted in Figure 10, for both the simple structured subdivision as well as the polygonal Voronoi-type one. In this figures, the subcell S c m under consideration would be colored red, while the subcells constituting N (S c m ) would be colored green. Let us emphasize that subcell S c m is also part of N (S c m ).

N (S c m ) =

Burgers equation with a smooth initial solution.

To highlight the efficiency of the developed APLSC-DG scheme in the non-linear case, let us first consider the 2D Burgers equation

(4.7a) (4.7b)    ∂ t u(x, t) + ∇ x F (u(x, t)) = 0, (x, t) ∈ [0, 1] 2 × [0, T ], u(x, 0) = u 0 (x), x ∈ [0, 1] 2 ,
where the convex flux function writes F (u) = 1 2 u 2 , u 2 t . As seen previously, starting from the smooth initial condition u 0 (x) = sin(2π (x, y)) on [0, 1] 2 , two stationary discontinuities form along the lines {(x, y) ∈ [0, 1] 2 , x + y = 1 2 } and {(x, y) ∈ [0, 1] 2 , x + y = 3 2 }. To emphasize how important a limiter or a correction oscillating the numerical is. Furthermore, the two shocks are absolutely not well resolved.

In Figure 29, the numerical solution obtained with the 6th-order APLSC-DG scheme is illustrated at time t = 0.5. We can see in this totally anisotropic triangles coarse grid, the corrected scheme quite accurately recovers the two straight line shocks, while ensuring a robust low oscillatory behavior. Now, to investigate once more if the cell subdivision type has any influence on the quality of the numerical solution obtained, we simulate this two-shocks Burgers test case with an even coarser grid with the four subdivisions depicted in Figures 14 and15 In the light of Figure 32, the four different subdivisions seem to produce equivalent results, which are further quite satisfactory considering the extremely coarse grid used. However, the use of the uniform structured cell partition, Figure 30(a), appears to capture in a sharper fashion the two straight line shocks. As we have shown that uniform subdivision, structured or Voronoi-type, lead to better results when APLSC-DG scheme is used, only those two subdivisions will be used in the remainder of the article.

KPP problem.

Before investigating the non-linear system case, we now turn our attention to non-linear conservation laws with non-convex fluxes. To this end, we consider the KPP problem proposed by Kurganov, Petrova, Popov (KPP) in [START_REF] Kurganov | Adaptive semi-discrete central-upwind schemes for non convex hyperbolic conservation laws[END_REF] to test the convergence properties of some WENO schemes in the context of non-convex fluxes. For this particular problem, we study the non-linear problem (4.7) where the flux function is given by F (u) = (sin(u), cos(u))

t . Considering the computational domain [-2, 2]×[-2.5, 1.5], the initial condition reads as follows

u 0 (x) =      7 π 2 if x < 1 2 , π 4 if x > 1 2 .
This test is very challenging to many high-order schemes as the solution has a two-dimensional composite wave structure, and as generally numerical methods fail to converge to the unique entropic exact solution.

In most cases, to be able to capture such rotation composite structure, very fine grids must be used. Here, by means of 6th-order uncorrected DG and then APLSC-DG scheme, we make use of an unstructured mesh made of 1054 triangular cells, which is very coarse in this quite complex situation. Results are displayed in Figure 33. Figure 33(a) illustrates the non-entropic behavior of DG schemes in this non-convex flux case, as well as the high oscillatory solution produced if no additional mechanism is used. The application of our APLSC technique permits to capture to correct entropic solution, while avoiding the apparition of spurious oscillations, as displayed in Figure 33(b). Furthermore, we can observe that although the coarseness of the grid used, the APLSC-DG scheme allowed to recover the two-dimensional vortex-like wave structure of the solution.

In Figure 33(c), subcells corrected during the different Runge-Kutta stages of the last time iteration are displayed, with different colors accordingly the amount of first-order correction applied. One can see how the NAD and PAD troubled subcell detection criteria accurately track the spiral discontinuity of the entropic exact solution. Here, a different colormap compared to Figure 6 has been used here for a better readability of the results. In Figure 34, we once more assess the impact of the subdivision on the quality of the results obtained the APLSC-DG method. As the uniform ones have proved to yield better results, we only compare here the uniform structured subdivision with uniform Voronoi-type one. Anew, the 4th-order scheme is used here. As one can see in Figure 34, similar results have been obtained regardless the type of cell subdivision. 4.3. 2D Euler system. To close this numerical application section, and assess once again the high capability of the a posteriori local subcell correction technique presented here, the non-linear system case will be now addressed. To this end, let us consider the 2D Euler compressible gas dynamics system of equations (4.8a)

(4.8b) (4.8c)                    ∂ t ρ + ∇ x q = 0, ∂ t q + ∇ x q ⊗ q ρ + p I d = 0, ∂ t E + ∇ x (E + p) q ρ = 0,
where the conserved variables ρ, q = ρ v and E respectively stand for the density, the momentum and the total energy, while v characterize the fluid velocity. The thermodynamic closure is given by the equation of state p = p(ρ, ε) where ε = E -1 2 ∥q∥ 2 ρ denotes the internal energy. In this paper, we make use of a gamma gas law, i.e. p = (γ -1) ε, where γ is the polytropic index of the gas.

Although the whole theory presented here has been introduced in the simple case of scalar conservation laws, the extension to the system case is perfectly straightforward. The only part which may differ is the troubled detection part. For the physical admissibility detection (PAD), we consider that a solution is physically admissible if the density and the internal energy are strictly positive. The use of other equations of state may lead a different convex set of admissibility, see [START_REF] Vilar | Positivity-preserving cell-centered Lagrangian schemes for multi-material compressible flows: From first-order to high-orders. Part I: The one-dimensional case[END_REF] for instance. For the numerical admissibility detection (NAD), the natural system counterpart would be to apply the previously introduced detection criteria to the Riemann invariants. However, in the non-linear system case, those quantities are not easy to get nor to manipulate. We could have use a linearized version of the Riemann invariants, as in [START_REF] Vilar | A discontinuous Galerkin discretization for solving the two-dimensional gas dynamics equations written under total Lagrangian formulation on general unstructured grids[END_REF] for instance, but for sake of simplicity we naively apply the NAD criterion to one of the conserved variable, similarly to what has been done in the 1D case [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. Here, we choose to work with the energy, as this physical quantity would be sensitive to any type of wave. Once again, the simple global Lax-Friedrichs numerical flux will be used in the remainder.

4.3.1. Sod shock tube problem.. We consider the extension of the classical Sod shock tube [START_REF] Sod | A survey of several finite difference methods for systems of non-linear hyperbolic conservation laws[END_REF] to the case of the cylindrical geometry. This problem consists of a cylindrical shock tube of unity radius. The interface is located at r = 0.5. At the initial time, the states on the left and on the right sides of the interface are constant. The left state is a high pressure fluid characterized by (ρ L 0 , p L 0 , v L 0 ) = (1, 1, 0), the right state is a low pressure fluid defined by (ρ R 0 , p R 0 , v R 0 ) = (0.125, 0.1, 0). The gamma gas law is defined by γ = 7 5 . The computational domain is defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, π 4 ]. We prescribe symmetry boundary conditions at the boundaries θ = 0 and θ = π 4 , and an outflow condition at r = 1. The exact solution consists of three circular waves, a shock followed by a contact discontinuity and rarefaction wave. The aim of this test case is then to assess the APLSC-DG scheme accuracy while ensuring a non-oscillatory behavior, and its ability to preserve the radial symmetry. In Figure 35, the 6th-order APLSC-DG scheme has been used on a very coarse anisotropic mesh made of 230 triangular cells. In the light of Figure 35(a), one can see how the radial wave structure has been accurately capture, even in this coarse mesh case. Figure 36, where the subcell mean values versus the subcell centroid radii x 2 + y 2 confirms this statement as the different points for a given radius do coincide. In Figure 35(b), subcells corrected during the different Runge-Kutta stages of the last time iteration are colored accordingly the amount of first-order FV correction used. It illustrates how this a posteriori correction procedure has been activated on zones corresponding to the solution loss of smoothness, meaning the left and right ends of the expansion fan, the contact discontinuity and the shock. One can also observe how the correction does operate locally inside the cell at a subcell scale, allowing the preservation of DG subcell high accurate resolution. Let us emphasize that this a posteriori correction procedure is not limited to the case of very high-order of accuracy on coarse grids. It also preforms very well at second or third order. See for instance Figure 37 where the third-order APLSC-DG numerical solution on 10571 triangular cells is depicted. 4.3.2. Sedov point blast problem.. We consider the Sedov problem for a point-blast in a uniform medium. An exact solution based on self-similarity arguments is available, see for instance [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF]. The initial conditions are characterized by (ρ 0 , p 0 , v 0 ) = (1, 10 -14 , 0), and the polytropic index is equal to 7 5 . We set an initial delta-function energy source at the origin prescribing the pressure in a control volume, yet to be defined, containing the origin as follows, p or = (γ -1) ε0 vor , where v or denotes the measure of the chosen control volume and ε 0 the total amount of release energy. By choosing ε 0 = 0.244816, as suggested in [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF], the solution consists of a diverging infinite strength shock wave whose front is located at radius r = 1 at t = 1, with a peak density reaching 6. The computational domain is defined in polar coordinates by (r, θ) ∈ [0, 1.2] × [0, π 4 ]. Similarly to the polar Sod shock tube problem, we prescribe symmetry boundary conditions at the boundaries θ = 0 and θ = π 4 , and an outflow condition at r = 1.2.

Regarding the control volume in which the delta-function energy will be dropped off, generally the cell containing the origin is considered. Here, to make this test case even more challenging, we choose to restrict the energy source only to the one subcell containing the origin. This means that initially, in one grid element the pressure in one subcell will be set to p or , while in the remainder of the cell the pressure will be 10 -14 . Let us further emphasize that generally in this test case, because one cannot simulate vacuum, the initial pressure is set to 10 -6 over the domain, except at the origin. Here, to make it once again more challenging, we set the initial pressure to 10 -14 .

We run this modified Sedov point blast problem with the sixth-order APLSC-DG scheme on a very coarse grid made of 271 triangular cells. In this particular case, the amount of total energy contained in the subcell located at the origin reaches 1947.5, while in the rest of the cell as well as in the remainder of the domain the total energy is set to 2.5E-14, see Figure 38. Any scheme lacking positivity-preserving property or a rigorous stabilization technique would fail solving this test problem. In Figure 38(b), one can see how the circular aspect of the shock has been accurately captured by the scheme, and the shock wave front is correctly located. This latter further goes different cells, enlightening the very precise subcell resolution of the APLSC-DG method. The numerical solution produced remains quite close to the one-dimensional self-similar exact solution, see Figure 39. Once again, this a posteriori local subcell correction procedure designed for DG schemes performs also very well for lower order methods, as depicted in Figure 40 where a finer grid made of 2894 cells has been used. Indeed, the numerical solution is very close to the one-dimensional analytical solution. In Figure 40(a), only the cell total energy means values are represented, and not the submean values as we generally do, for a better readability of the results in this finer grid context.

4.3.3.

The forward-facing step problem.. We now consider the forward facing step problem, which has been initially introduced by A. Emery in [START_REF] Emery | An evaluation of several differencing methods for inviscid flow problems[END_REF], and further studied by P. Woodward and P. Colella in [START_REF] Woodward | The numerical-simulation of two-dimensional fluid-flow with strong shocks[END_REF]. This challenging test case consists in a Mach 3 flow in a 3 units in length and 1 unit in width wind tunnel. Initially, the tunnel is filled with a gamma gas law with γ = 7 5 , which everywhere has density ρ 0 = 1.4, pressure p 0 = 1 and velocity v 0 = (3, 0)

t . The 0.2 high step being located at x = 0.6, the computational domain is then ([0, 3] × [0, 1]) \ ([0.6, 3] × [0.2, 1]). Gas with this density, pressure and velocity is continually fed in from the left-hand boundary. Let us emphasize that unlike as it is generally done, we did not refine the mesh near the corner, see Figure 41, nor modify in any way our APLSC-DG scheme. In Figure 41, the numerical solution obtained by means of 6th-order APLSC-DG scheme on an unstructured grid made of 680 cells is displayed. Let us note that despite the coarseness of the mesh used, the shocks and the rarefaction fan created around the corner are quite well resolved, while ensuring an low oscillatory behavior. This result demonstrates once again the high capability of the presented a posteriori local subcell correction combined with high-order discontinuous Galerkin methods.

Conclusion.

The paper aims at presenting the two-dimensional unstructured extension of the new correction technique of DG schemes introduced in [START_REF] Vilar | A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. This a posteriori procedure relies on the expression of DG methods as a FV-like scheme on a subgrid. By means of this theoretical part, we modify at the subcell level the so-called reconstructed fluxes only where the uncorrected DG scheme has failed. Consequently, only very few subcells require this particular treatment. In this paper, a new version of the correction procedure is also introduced, where a convex blending of high-order DG reconstructed fluxes and first-order FV fluxes is applied in the vicinity of troubled zones. For the remaining subcells, the submean values obtained through the uncorrected DG method are kept, as they have been detected as admissible by troubled zone criteria. This correction procedure allows us to retain the very precise subcell resolution of DG schemes, along with addressing the issues of spurious oscillations or non-entropic behavior. A wide number of test cases on different problems have been used to depict the good performance and robustness of the presented correction technique. Different types of cell subdivision have been compared.

In the future, we intend to extend this a posteriori correction technique to moving grid configurations, both ALE and Lagrangian formalisms, as well as the case of curvilinear meshes. We also plan to adapt this local subcell reconstructed flux correction framework to the a priori paradigm, by means of the Flux-Corrected Transport (FCT) methodology, in order to obtain an automatic very high-order and property preserving scheme.

Fig. 1 . 1 S 2 S 3 ( 8 (c) P 3 .Fig. 2 . 3 .

 11238323 Fig. 1. Examples of subdivision for a triangular cell and a P 3 DG scheme.

  [START_REF] Gassner | A Skew-Symmetric Discontinuous Galerkin Spectral Element Discretization and Its Relation to SBP-SAT Finite Difference Methods[END_REF] whereA c ∈ M N k ×N c f , defined as (A c ) mp = ε c mp ,stands for the adjacency matrix, and where D c = diag |S c 1 |, . . . , |S c N k | is the subcells volume matrix and (B c ) m = ∂S c m ∩∂ωc F n dS carries the cell boundary contribution. Finally, introducing the matrix Q c = D c P c and substituting equation (2.13) into (2.20), one gets

  ) where v min x\y, p and v max x\y, p are respectively defined as v min x\y, p = min v∈Vp v c x\y (x p ) and v max x\y, p = max v∈Vp v c x\y (x p ). (3.6)

  Structured subdivision. (b) Voronoi-type subdivision.

Fig. 3 .

 3 Fig. 3. Original correction of the DG reconstructed flux.

1 .. 3 .

 13 Compute the candidate solution u c, n+1 h by means of the uncorrected DG scheme. 2. Project u c, n+1 h through (3.2) to get the N k submean values u c, n+1 m Check u c, n+1 m through the troubled zone detection criteria plus relaxation . 4. If all submean values u c, n+1 m are admissible, go further in time. Otherwise, modify the corresponding reconstructed flux value on the troubled subcells faces through a first-order numerical flux as following

6 .. 7 .

 67 By means of P -1 c and equation (3.2), get the new corrected polynomial solution u c, n+1 h Return to point 3.

  a) Structured subdivision. (b) Voronoi-type subdivision.

Fig. 4 .

 4 Fig. 4. New correction of the DG reconstructed flux.

  t = 0.5.

Fig. 5 .

 5 Fig. 5. Subcell mean values of the entropic weak solution of the Burgers equation with u 0 (x) = sin(2π (x + y)).

Fig. 6 .

 6 Fig. 6. Comparison between original and new correction procedure: corrected subcells.

Fig. 7 .

 7 Fig. 7. Comparison between original and new correction procedure: corrected subcells -zoom on [0.65, 0.9] 2 .

Fig. 8 .Fig. 9 .

 89 Fig. 8. Comparison between original and new correction procedure: subcell mean values.

(

  ||F ′ (u c m )|| 2 ), and where the cell and subcell characteristic lengths d c and d c m are defined as follows d c = |ω c |

  Voronoi-type subdivision.

Fig. 10 .

 10 Fig. 10. Neighboring subcells set N (S c m ) for the NAD criterion in the linear case: subcell S c m is colored red while the subcells in N (S c m ) are colored green.

Fig. 11 .

 11 Fig.11. Linear advection a smooth signal with a 6th corrected DG scheme, on a 5x5x4 grid after 1 period.

  APLSC-DG: min=8.5E-6, max=1.

Fig. 12 .

 12 Fig. 12. 6th-order DG solutions for the linear advection of a crenel signal on 576 cells after one period.

20 -Fig. 13 .

 2013 Fig. 13. 6th-order solutions for the crenel advection case on 576 cells: submean values versus (x + y -1) coordinate.

  (a) Equidistant boundary points. (b) Gauss-Lobatto boundary points.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Examples of structured subdivisions for a triangular cell and a P 3 DG scheme.

  , the subdivision is obtained as follows: first, we define a triangulation of the element by joining the P k Lagrangian nodes. Then, the centroid of each triangle and the midpoint of the edges form a set of N k polygonal subcells. In Figure15(a), the previous procedure is modified to yield a more uniform subdivision by setting equidistant cell boundary points. Now, to make sure that cell subdivision does not have any impact on the numerical solution, we display in Figures 16, 17 and 18 the solutions obtained through subcell finite volume scheme with high-order reconstructed fluxes as numerical fluxes, equation (2.40), with the four different cell subdivisions previously introduced. In Figure18, where submean values versus (x + y -1) coordinate are displayed, it is clear that Non-uniform structured subdivision.

Fig. 16 .

 16 Fig. 16. 4th-order DG solutions for the crenel signal advection on 576 cells after five periods: structured subdivision.

  Non-uniform polygonal subdivision.

Fig. 17 .Fig. 18 .

 1718 Fig. 17. 4th-order DG solutions for the crenel signal advection on 576 cells after five periods: polygonal subdivision.

  Non-uniform structured subdivision.

Fig. 19 .

 19 Fig. 19. 4th-order APLSC-DG solutions for crenel advection on 576 cells after five periods: structured subdivision.

Fig. 20 .

 20 Fig. 20. 4th-order APLSC-DG solutions for crenel advection on 576 cells after five periods: polygonal subdivision.

Fig. 21 .

 21 Fig. 21. 4th-order APLSC-DG solutions for the crenel advection case on 576 cells using different cell subdivisions: submean values versus (x + y -1) coordinate.

  t = 2 π: min=6.7E-10, max=0.97.

Fig. 22 .

 22 Fig. 22. 6th-order APLSC-DG solution for the rigid rotation case on 576 cells.

  APLSC-DG (b) Solution profile for x = 0.25.

Fig. 23 .

 23 Fig. 23. 6th-order APLSC-DG solution for rigid rotation on 576 cells after one full rotation: solution profiles.

Fig. 24 .

 24 Fig. 24. 4th-order APLSC-DG solutions for rigid rotation on 576 cells after five full rotations: structured subdivision.

  Non-uniform polygonal subdivision.

Fig. 25 .

 25 Fig. 25. 4th-order APLSC-DG solutions for rigid rotation on 576 cells after five full rotations: polygonal subdivision.

Fig. 26 .

 26 Fig. 26. 4th-order APLSC-DG solutions for rigid rotation on 576 cells after five full rotations: solution profiles.

  Voronoi-type subdivision.

Fig. 27 .

 27 Fig. 27. Neighboring subcells set N (S c m ) for the NAD criterion in the non-linear case: subcell S c m is colored red while the subcells in N (S c m ) are colored green.

  Submean values versus (x + y -1) coordinate.

Fig. 28 .

 28 Fig. 28. 6th-order uncorrected DG solution for 2D Burgers equation on a 576 cells mesh at t = 0.5.

  .

  Submean values versus (x + y -1) coordinate.

Fig. 29 .

 29 Fig. 29. 6th-order APLSC-DG solution for 2D Burgers equation on a 576 cells mesh at t = 0.5.

Fig. 30 .

 30 Fig. 30. 4th-order APLSC-DG solutions for 2D Burgers equation on 242 cells at t = 0.5: structured subdivision.

  Non-uniform polygonal subdivision.

Fig. 31 .Fig. 32 .

 3132 Fig. 31. 4th-order APLSC-DG solutions for 2D Burgers equation on 242 cells at t = 0.5: polygonal subdivision.

  DG solution map.

  APLSC-DG solution map.

  Corrected subcells map.

Fig. 33 .

 33 Fig. 33. 6th-order APLSC-DG solution for the KPP problem on a 1054 cells mesh at time t = 1.

  Voronoi-type subdivision.

Fig. 34 .

 34 Fig. 34. 4th-order APLSC-DG solution for the KPP problem on a 1054 cells mesh at time t = 1: subdivisions comparison.

Fig. 35 .Fig. 36 .

 3536 Fig. 35. 6th-order APLSC-DG solution for the polar Sod shock tube problem on 230 cells.

  APLSC-DG (b) Submean values versus radii.

Fig. 37 .

 37 Fig. 37. 3rd-order APLSC-DG solution for the polar Sod shock tube problem on 10571 cells.

Figure 37

 37 Figure37depicts how the numerical solution produced by our third-order APLSC-DG scheme respects the radial symmetry of the flow, and tends to converge to the exact solution solution in this quite refined mesh context.

Fig. 38 .Fig. 39 .

 3839 Fig. 38. 6th-order APLSC-DG solution for the Sedov problem on 271 cells: subcell mean total energy values.

  Subcell mean density versus radii.

Fig. 40 .Fig. 41 .

 4041 Fig. 40. 3rd-order APLSC-DG solution for the Sedov point blast problem at t = 1 on 2894 cells.

Table

	1 40	4.03E-11	-	4.52E-11	-	1.50E-10	-

Table 2

 2 Projection matrix condition number for different orders and subdivisions.

	subdiv.	1	4	10.91 31.75
	Non-unif. struct. subdiv.	1	4	9.52 29.28
	Unif. polyg. subdiv.	1 2.87 8.73 27.89
	Non-unif. polyg. subdiv.	1 2.87 8.19 26.94

Acknowledgment. R.A. was partially funded by SNF project "Structure preserving and fast methods for hyperbolic systems of conservation laws" number 200020 204917.