A posteriori subcell finite volume limiter for general PNPM schemes: applications from gasdynamics to relativistic magnetohydrodynamics

In this work, we consider the general family of the so called ADER PNPM schemes for the numerical solution of hyperbolic partial differential equations with \textit{arbitrary} high order of accuracy in space and time. The family of one-step PNPM schemes was introduced in [Dumbser et al., JCP, 2008] and represents a unified framework for classical high order Finite Volume (FV) schemes (N=0), the usual Discontinuous Galerkin (DG) methods (N=M), as well as a new class of intermediate hybrid schemes for which a reconstruction operator of degree M is applied over piecewise polynomial data of degree N with M>N. In all cases with M>= N>0 the PNPM schemes are linear in the sense of Godunov, thus when considering phenomena characterized by discontinuities, spurious oscillations may appear and even destroy the simulation. Therefore, in this paper we present a new simple, robust and accurate a posteriori subcell finite volume limiting strategy that is valid for the entire class of PNPM schemes. The subcell FV limiter is activated only where it is needed, i.e. in the neighborhood of shocks or other discontinuities, and is able to maintain the resolution of the underlying high order PNPM schemes, due to the use of a rather fine subgrid of 2N+1 subcells per space dimension. The paper contains a wide set of test cases for different hyperbolic PDE systems, solved on adaptive Cartesian meshes (AMR) that show the capabilities of the proposed method both on smooth and discontinuous problems, as well as the broad range of its applicability. The tests range from compressible gasdynamics over classical MHD to relativistic magnetohydrodynamics.


Introduction
In this work we want to improve the family of high order accurate ADER P N P M schemes first introduced in [39,36] for the solution of hyperbolic partial differential equations. In this family of schemes the discrete solution is represented in space through high order piecewise polynomials of degree N at each timestep; the data are then evolved in time through a space-time reconstruction procedure of order M . The reconstruction procedure is divided into two steps: concerning the spatial reconstruction, we employ a classical WENO reconstruction in the case of pure finite volume schemes (N = 0), a reconstruction procedure based on L 2 projection that is linear in the sense of Godunov for N > 0 and M > N , and in the case or pure DG schemes (N = M ) the reconstruction reduces to the identity operator; concerning the reconstruction in time, we employ a novel variant of the ADER approach of Toro and Titarev, see [105,109,106,110,26], based on an element-local spacetime Galerkin predictor, see [39]. In practice, we can see the Finite Volume (FV) schemes of order M as a particular case of P N P M methods when N = 0, and also the Discontinous Galerkin (DG) methods are included in this family when choosing N = M .
Furthermore, this family contains another important class of hybrid or reconstructed DG schemes when taking N > 0, M > N , which are the main object of study of this paper. Indeed, they offer many advantages, in particular their good compromise between cost and resolution. In fact, data are represented with polynomials of order N , so more accurately with respect to FV methods, but without the expensive cost of a full DG representation of approximation degree M ; also the CFL stability constraint that limits the timestep size of any explicit scheme, only depends on N and not on M , allowing for larger timesteps once the desired order of accuracy has been fixed, see [39]. Last but not least, for N > 0 the P N P M schemes require a much smaller reconstruction stencil than comparable finite volume schemes of degree M . The nominal order of accuracy of the scheme is given by M + 1 so it can be at least in principle arbitrary high.
The family of reconstructed DG schemes, which is similar to the P N P M framework, was forwarded independently by Luo et al. in a series of papers, see e.g. [84,85,65,113,115,112] and references therein. At this point we also highlight that the use of reconstruction and filtering operators as a post-processor for improving the accuracy of DG schemes goes back to work of Ryan et al., see [96,69,95,74,87]. Other related work on reconstruction-based DG schemes can be found in [79,28,114,17].
Moreover, the ADER P N P M family provides a useful framework for code developers because it allows to include in a unique code both types of standard discretization methods for hyperbolic PDE (FV and DG schemes), together with the new class of intermediate hybrid schemes for M > N > 0. It is then possible to let it up to the user to decide whether for a particular application the use of a robust finite volume approach (N = 0), a very accurate DG scheme (N = M ), or a less expensive but still very accurate intermediate P N P M method with M > N > 0, is the most appropriate.
The main drawback so far of the intermediate P N P M schemes with M > N > 0, as presented in [39], is that they are linear in the sense of Godunov [61], hence not well suited for dealing with discontinuous problems. For this reason, here we propose a new simple, robust and accurate limiting strategy that is able to stabilize the entire class of P N P M schemes in such a way that they can be employed for the numerical solution of hyperbolic equations with discontinuous solutions, which may arise even when starting with smooth initial conditions. Moreover, the new limiter does not substantially deteriorate the benefits of P N P M schemes in terms of computational cost and accuracy of the original unlimited schemes. To the best knowledge of the authors, this is the first time that an a posteriori subcell finite volume limiter is proposed for general P N P M schemes with M > N > 0. So far, only the cases N = 0 and M = N were covered in [39] and [50], respectively.
Our limiter is based on the MOOD approach [29,34,35], which has already been successfully applied in the framework of ADER finite volume schemes [82,22,21] and Discontinous Galerkin finite element schemes, see [50,120,45,25,72]. Specifically, the numerical solution is checked a posteriori for nonphysical values and spurious oscillations, and if it does not satisfy all admissibility detection criteria, given by both physical and numerical requirements, in a certain cell, that cell is marked as troubled. Then, instead of applying a limiter to the already computed solution, the solution is locally recomputed with a more robust scheme in the troubled cells, relying either on a second order TVD scheme, as proposed for pure DG schemes in [45,19,102], or on a higher order ADER-WENO finite volume method as employed in [50,120,23,92,94]. Moreover, this second computation is performed on a finer subgrid generated within each troubled cell; the subcell approach is employed in order to maintain the high resolution of the initial P N P M scheme even when passing to a less accurate (but more robust) FV scheme. For the given reasons our limiter is called a posteriori subcell finite volume limiter.
Finally, for a complete review of ADER P N P M schemes we refer to the recent paper [25], where a complete introduction traces the historical developments of these methods up to its latest evolutions.
The rest of the paper is organized as follows. After an introduction of the class of physical phenomena that can be discretized with the proposed numerical method and the structure of our data representation, we present the family of ADER P N P M schemes in Section 2. In particular, we describe the reconstruction procedure in space, see Section 2.3, and in time see Section 2.4; these procedures provide a high order reconstructed polynomial of degree M in space and time that will be used in the final one-step update formula given in Section 2.5. Then, Section 2.6 is dedicated to our a posteriori subcell FV limiter, which in addition can be combined with mesh adaptation techniques as described in Section 2.7.
Next, in Section 3 we present a large set of numerical results that shows the order of convergence of our scheme for smooth solutions and their capability of dealing with discontinuities, i.e. their robustness and resolution. We also compare the hybrid reconstructed schemes with pure DG schemes in order to show the resulting gain in terms of computational cost. Finally, we close the paper with some remarks and an outlook to future works in Section 4.

Numerical method
In this Section we carefully describe the a posteriori subcell finite volume limiter for general P N P M schemes, showing its simplicity, accuracy, robustness and versatility thanks to the following key ingredients: the use of the unified P N P M framework for finite volume (FV), discontinous Galerkin (DG) and hybrid reconstructed DG schemes allows the user to decide freely which combination of N and M is the better choice for a particular application; the ADER space-time predictor-corrector formalism allows the implementation of a truly arbitrary high order accurate fully discrete one-step scheme that needs only one MPI communication per time step within a parallel HPC implementation, see Section 2.4; the a posteriori subcell finite volume limiter avoids spurious oscillations of high order P N P M schemes without affecting the resolution of the underlying method, see Section 2.6; the adaptive mesh refinement (AMR) technique allows to use a fine grid only where necessary, resorting to cheaper coarse grids in smooth regions of the solution, see Section 2.7.

Governing PDE system
We consider a very general formulation of the governing equations in order to model a wide class of physical phenomena, namely all those which are described by hyperbolic systems of conservation laws that can be cast into the following form, where x = (x, y, z) is the spatial position vector, d is the number of space dimensions, t represents the time, Q = (q 1 , q 2 , . . . , q m ) T is the vector of conserved variables defined in the space of the admissible states Ω Q ⊂ R m and is the non-linear flux tensor. This kind of system (1) is said to be hyperbolic if for all directions n = 0 the matrix A n = ∂F/∂Q · n has m real eigenvalues and a full set of m linearly independent eigenvectors. Examples of hyperbolic equations are the Euler equations of gasdynamics, the Shallow Water equations [27,104] and many multiphase models [4,37,59] used in fluid mechanics, the magnetohydrodynamics system (MHD) for plasma flow [9,6], the unified first order hyperbolic formulation of continuum mechanics by Godunov, Peshkov and Romenski (GPR) [63,90,62,46,47,41] as well as the special and general relativistic formulations of MHD, see e.g. [5,118,10,3,32,53], or for the Einstein field equations (CCZ4) [1,2,40,42]. We will test the method proposed in this paper on some of those systems in order to verify its applicability in different physical domains.
2.2 Domain discretization and high order data representation (order N ) On grid level = 0 the computational domain Ω is discretized with a uniform Cartesian grid, called main grid or the level zero grid, For each element we define a reference frame of coordinates ξ = (ξ, η, ζ) linked to the Cartesian coordinates x = (x, y, z) of Ω ijk by (3) Then, we represent the conserved variables Q of (1) in each cell Ω i by a d−dimensional tensor product of piecewise polynomials of degree N where ϕ (ξ) are nodal spatial basis functions given by the tensor product of a set of Lagrange interpolation polynomials of maximum degree N such that where ξ m GL are the set of (N + 1) d Gauss-Legendre (GL) quadrature points obtained by the tensor product of the GL quadrature points ξ m GL , η m GL , ζ m GL in the unit interval [0, 1], see [103].
The discontinuous finite element data representation in (4) leads naturally to i) a Discontinuous Galerkin (DG) scheme if N > 0 and N = M , where the desired order of accuracy M already coincides with the degree N of the polynomial approximating the data (M = N ), so that high order of accuracy in space can be obtained without the use of any spatial reconstruction operator, and to ii) a Finite Volume (FV) scheme in the case N = 0. This indeed means that for N = 0 we have ϕ (ξ) = 1 with = 0, and (4) reduces to the classical piecewise constant data representation that is typical of finite volume schemes, where the only degree of freedom per element is the usual cell averagê u 0 . In this case the order of accuracy M in space will be obtained through the reconstruction procedure described in next Section 2.3. However, iii) also a family of hybrid reconstructed Discontinuous Galerkin methods is included in this representation, where a Hermite-type reconstruction of degree M > N is performed on cell data represented by polynomials of degree N , see the next Section 2.3.
Thus, within the general P N P M formalism one can simultaneously deal with arbitrary high order FV and DG schemes and reconstructed hybrid methods inside a unified framework, with only very few differences between the different schemes (substantially the reconstruction procedure and the type of limiter).

High order spatial reconstruction (order M )
In the framework of P N P M schemes, M indicates the highest polynomial approximation degree used for the representation of the discrete solution within the method. Hence, in this Section we describe the reconstruction procedure that is needed to obtain approximation degree M in space from an underlying data representation u h (x, t n ) of lower or equal degree N ≤ M , i.e. the procedure that generates a spatially high order accurate reconstruction polynomial where we formally employ the same nodal basis functions for the reconstruction and for the data representation, see (4). However, note that when M = N of course ψ l (x, t n ) does not coincide with ϕ l (x, t n ), since the polynomial degree and the positions of the GL points are not the same.
For the sake of a uniform notation, when M = N , we trivially impose that the reconstruction polynomial is given by the DG polynomial, i.e. w h (x, t n ) = u h (x, t n ), which automatically implies that in the case N = M the reconstruction operator is simply the identity.
In the other cases, we employ a polynomial reconstruction procedure implemented in a dimension by dimension fashion in order to compute the coefficientsŵ ,i in (6). To better follow the following reasoning we refer the reader Since we are employing nodal basis functions, we can represent the available information at each stage of our P N P M scheme in each cell by a symbol located at a certain GL point inside the cell. In a P 1 P 2 scheme u h is represented by a P N = P 1 polynomial, so we have (N + 1) 2 = 4 information in each cell (a, the green circles). By selecting the N + 1 = 2 information along the same horizontal section in Ω ij and in its two immediate neighbors Ω i−1,j , Ω i+1,j (b), we have enough information (3(N + 1) = 6 > 3 = M + 1) in order to reconstruct a P M = P 2 polynomial in x-direction (c); then we have to repeat the same procedure for each N + 1 = 2 horizontal section of cell Ω ij (d-e). In this way we obtain our reconstructing polynomial in the x-direction, represented by (M + 1)(N + 1) = 6 information (f, the blue crosses).
also to the Figures 1, 2, and 3. Focusing on the reconstruction procedure along the x-direction, given an element Ω i = Ω ijk , we write the reconstruction polynomial in x-direction w x h in terms of one dimensional basis functions as Then, we integrate on a set S x of neighbors of Ω i in x-direction, obtaining an algebraic system for the polynomial coefficientsŵ 1 ,r2,r3,i (one for each  After having performed the reconstruction along the x direction, in each cell we have (M + 1)(N + 1) = 6 information (a, the blue crosses). Now, by selecting (N + 1) = 2 information along the same vertical section in Ω ij and in its two immediate neighbors Ω i,j−1 , Ω i,j+1 (c), we have enough information (3(N + 1) = 6 > 3 = M + 1) in order to reconstruct a P M = P 2 polynomial in y-direction (d); then we have to repeat the same procedure for each M + 1 = 3 vertical section of cell Ω ij (e-h). In this way we obtain our final P M = P 2 reconstructing polynomial for the cell Ω ij (i, the red crosses).  The available data (a, green circles) are provided by the P N = P 2 polynomial u h . By selecting (N + 1) = 3 information along the same horizontal section in Ω ij and in its two immediate neighbors Ω i−1,j , Ω i+1,j (b), we have enough information (3(N +1) = 9 > 5 = M +1) in order to reconstruct a P M = P 4 polynomial in xdirection (c); then the same procedure has to be repeated for each N + 1 = 3 horizontal section of cell Ω ij obtaining (d), and finally for each cell of the domain. At this point in each cell we have (M + 1)(N + 1) = 15 information (d, the blue crosses) and by selecting (N + 1) = 3 information along the same vertical section in Ω ij and in its two immediate neighbors Ω i,j−1 , Ω i,j+1 (e), we have enough information to reconstruct a P M = P 4 polynomial in ydirection (f); then we have to repeat the same procedure for each M + 1 = 5 vertical section of cell Ω ij (g-h). In this way we obtain our final P M = P 5 reconstructing polynomial for the cell Ω ij (i, the red crosses).
where S x contains the neighbors along the x axis, i.e Ω m = Ω mjk such that m ∈ {i − 1, i, i + 1}. Note that we are using a stencil made always of only 3 elements in each direction, thus a very compact one. Indeed, a stencil composed of 3 elements is enough for any P N P M scheme with M ≤ 3N + 2, because any P N cell contains N + 1 degrees of freedom, thus 3 cells provide 3N + 3 degrees of freedom, which are sufficient for a polynomial reconstruction of degree up to 3N + 2. Moreover, when the provided information are more than the minimum required, the system (9) results to be overdetermined; so, to solve it, we employ a constrained least-squares technique (CLSQ) [44], i.e. we impose that the reconstructed polynomial satisfies exactly. In other words, all moments of the reconstructed solution w h and the original solution u h up to degree N must coincide exactly within cell Ω i and match on the remaining stencil elements in the least-square sense.
To complete the reconstruction polynomial, we now repeat the above procedure in the y-direction, so we write the reconstruction polynomial in terms of one-dimensional basis functions as and we solve the algebraic system 1 ∆y n with S y being the set of the neighbors along the y axis, i.e Ω n = Ω ink such that n ∈ {j − 1, j, j + 1}. Again, for overdetermined systems we impose that the reconstruction exactly satisfies And then the same procedure can be repeated along the z axis by looking for the unknown coefficientsŵ 1, 2 , 3 ,i of and solving the algebraic system S z being the set of the neighbors along the z axis, i.e Ω p = Ω ijp such that p ∈ [k − 1, k, k + 1]. For overdetermined systems, the constraint reads Finally, the coefficientsŵ 1 , 2, 3,i represent theŵ ,i of (6) that give us the desired polynomial representation of order M in space.
We would like to emphasize that the reconstructed P N P M schemes with N > 0 are very compact because for the reconstruction they need a much smaller stencils than classical finite volume schemes and that for regular Cartesian meshes, the coefficients of the above constraint least squares systems depend only on the choice of the basis functions, hence the integrals can be precomputed once and for all on the reference element before starting the simulation.
Furthermore, in the specific case N = 0, i.e. when the P N P M reduces to a FV scheme, the above polynomial reconstruction procedure must be made nonlinear; this can be easily done, for example, by adopting the WENO strategy specifically described in the context of P N P M type schemes (thus with the same notation adopted here) on Cartesian meshes in [49,25]. We recall that the nonlinearity introduced through ENO/WENO type procedures essentially avoids the spurious oscillations typical of high order linear schemes modeling discontinuous processes see [61], and was already introduced in the 80s and subsequently largely developed [67,66,99,70,8,121,101]. Due to the already exhaustive literature available on FV schemes, here, for what concerns the strategies that guarantee robustness on discontinuities, we focus on P N P M schemes only with N > 0: indeed, it is for those schemes that we propose in this work a new strategy, i.e. the new a posteriori subcell FV limiter described in Section 2.6.

High order in time via a local space-time Galerkin predictor
We recall that high order of accuracy in space is provided by the piecewise polynomial data representation w h of (6), obtained in the previous Section 2.3. Now, in order to achieve also high order of accuracy in time, relying on the ADER predictor-corrector approach, we need to compute the so-called space-time Galerkin predictor, i.e. a space-time polynomial q h of degree M in (d + 1)-dimensions (d for the space plus 1 for the time) which takes the following form where again θ (ξ, τ ) is given by the tensor product of Lagrange interpolation polynomials ϕ (ξ(x)) ϕ τ (τ ), with ξ(x) given by (3) and the mapping for the time coordinate given by t = t n +τ ∆t, τ ∈ [0, 1]. This high order polynomial in space and time will serve as a predictor solution, only valid inside Ω i ×[t n , t n+1 ], to be used for evaluating the numerical fluxes and the sources when integrating the PDE in the final corrector step of the ADER scheme, see Section 2.5.
In order to determine the unknown coefficientsq of (16) we search q h such that it satisfies a weak form of the governing PDE (1) integrated in space and time locally inside each where the first term is integrated in time by parts exploiting the causality principle (upwinding in time) and w h (x, t n ) is the known initial condition at time t n . Now, the system (18), which contains only volume integrals to be calculated inside Ω i and no surface integrals, can be solved via a simple discrete Picard iteration for each element Ω i , and there is no need of any communication with neighbor elements. Indeed, the so-called predictor step consists in a local solution of the governing PDE (1) in the small, see [66], inside each spacetime element Ω i × [t n , t n+1 ]. It is called local because it is obtained by only considering cell Ω i with initial data w h , the governing equations (1) and the geometry, without taking into account any interaction between Ω i and its neighbors. We also want to emphasize that this procedure is exactly the same whatever N and M are.
We recall that this procedure has been introduced for the first time in [39] for unstructured meshes, it has been extended for example to moving meshes in [18] and to degenerate space time elements in [57]; finally, its convergence has been formally proved in [25].
2.5 High order fully-discrete one-step ADER P N P M scheme Last, the update formula of our ADER P N P M scheme is recovered starting from the weak formulation of the governing equations (1) (where the test functions ϕ k coincide with the basis functions ϕ of (5)) we then substitute Q with (4) at time t = t n (the known initial condition) and at t = t n+1 (to represent the unknown evolved conserved variables), and with the high order predictor q h previously computed for t ∈ [t n , t n+1 ], obtaining The use of q h allows to compute the integrals appearing in (20) with high order of accuracy in both space and time.
The boundary fluxes F ·n are obtained by a Riemann solver, thus providing the coupling between neighbors, which was neglected in the predictor step. In particular, in this work we will employ three types of standard fluxes, namely the Rusanov flux and the HLL flux, whose description can be found in [108], and the HLLEM flux for which we refer to [51,38]. For the sake of completness, we report here the expression of the Rusanov flux that reads as follows where s max is the maximum eigenvalue of the system matrices A(q + h ) and We remark also that due to the discontinuous character of q h at the interfaces ∂Ω i , F ·n is computed through a numerical flux function evaluated over the boundary-extrapolated data q − h and q + h (i.e the predictors q h of two neighbors elements evaluated at the common interface).
Finally, we stress again that the update procedure in (20) is the same whatever N and M are, and allows the contemporary evolution of all the (N + 1) d degrees of freedom of u h .

CFL stability constraint
A very important feature of P N P M schemes is linked to the CFL stability constraint. Since this family of scheme is explicit, the time step ∆t has to be computed according to a (global) Courant-Friedrichs-Levy (CFL) stability condition given by where h min is the minimum characteristic mesh-size, |λ max | is the spectral radius of the system matrix A and the maximum admissible CFL P N P M number is given in Table 1. In the above formula we wanted also to recall the classical CFL condition of Runge-Kutta DG schemes (the one written on the right, with CFL < 1) which is just a bit less restrictive than the one needed for ADER P N P M schemes, but easier to remember and helpful in justifying the stability of our subcell limiter, see formula (29). Furthermore, we would like to emphasize that it is the degree N of the data representation that governs the stability of the method and not the polynomial degree M of the reconstruction operator. Hence, the reconstructed hybrid methods with N > 0 and M > N allow for larger time steps than the pure DG methods (N = M ) of the same order always maintaining a superior resolution with respect to FV schemes (N = 0), fact that further justifies the interest in their development.

A posteriori subcell finite volume limiter
Up to now, the presented P N P M scheme is high order accurate in space and time and, formally, the differences between the FV case (N = 0) the pure DG case (N = M ) and the hybrid reconstructed case (N > 0, M > N ) are basically only due to the procedure for achieving high order of accuracy in space, which is obtained through a WENO reconstruction in the FV case, a linear reconstruction in the hybrid case and is automatic by construction for DG, see Section 2.3. But this is actually a major difference, because the WENO operator provides a non-linear stabilization of the FV scheme, while the P N P M schemes with N > 0 presented so far are unlimited and, as such, they are affected by the so-called Gibbs phenomenon, i.e. oscillations are likely to appear in presence of shock waves or other discontinuities. These oscillations can be explained by the Godunov theorem [61], because in this case the scheme is linear in the sense of Godunov.
As a consequence, a limiting technique is required. Our strategy is described in detail below and it will be applied whenever N > 0.
First, we need to consider the numerical solution computed so far u n+1 h only as a candidate solution: we denote it with u n+1, * h (x, t n+1 ). Then, following [50,120,119,54,25,72], each element Ω i is divided into N ω = (2N + 1) d equal non-overlapping subgrid cells ω i,α whose volume is denoted by |ω i,α |; for any cell we define the corresponding subcell average of the and the candidate subcell averages where P(u h ) is the L 2 projection operator into the space of piecewise constant cell averages. Now, we have to mark the troubled cells, i.e. we have to identify those cells where the solution found through the P N P M scheme cannot be accepted because it may lead to spurious oscillations. Thus, the candidate solution v n+1, * h is checked against a set of detection criteria. Here we follow the criteria described in [19], however also other specific physical bounds or more elaborate choices as those of invariant domain preserving methods [64] could be considered.
First, we require that the computed solution is physically acceptable, i.e. that it belongs to the phase space of the conservation law being solved. For instance, if the compressible Euler equations for gas dynamics are considered, density and pressure should be positive and in practice we require that they are greater than a prescribed tolerance ε = 10 −12 . Then, the solution should verify a relaxed discrete maximum principle (DMP) where V(Ω i ) is the set containing all the neighbors of Ω i sharing a common node with Ω i , and δ is a parameter which, according to [19,50,120], reads with δ 0 = 10 −5 and ε = 10 −4 . If a cell does not fulfill the detection criteria in all its subcells, then it is marked as troubled. It is possible that some false positive activations of the limiter occur; however these local effects do not reduce the overall quality of the simulation thanks to the highly accurate limiter procedure adopted on troubled cells.
Then only on these troubled cells we apply either a second-order accurate MUSCL-Hancock TVD finite volume scheme with minmod slope limiter [108] (in particular in presence of strong shock waves or low density atmospheres), or a more accurate ADER-WENO FV scheme [50,49] that better captures local extrema. In this way we can re-compute the solution in order to evolve the cell averages v n i,α in time and obtain v n+1 i,α . Note that, due to the fact of applying a high order scheme and to do so on a subgrid instead that on the main grid, the subcell average representation given by v n+1 i,α maintains the high resolution of the underling P N P M scheme. Indeed now, we can recover from these cell averages a polynomial u n+1 h of degree N ; this is done by applying a reconstruction operator R such that which is conservative on the main cell Ω i thanks to the additional linear constraint Moreover, the projection operator P in (23) and the reconstruction operator R in (27) satisfy the property P · R = I, with I being the identity operator. However, we have to remark that the reconstruction operator (27)-(28) might still lead to an oscillatory solution, since it is based on a linear unlimited least squares technique. If this is the case, the cell Ω i will be marked again automatically as troubled during the next timestep t n+2 , therefore the same finite volume subcell limiter will be used again in that cell and in particular the subcell averages from which to start as initial data at time t n+1 will be the v n+1 i,α kept in memory from the previous limited step. Furthermore, if a cell Ω i is acceptable but has at least one troubled neighbor in V(Ω i ), then we cannot accept its candidate solution u n+1, * h (x, t n+1 ) as it is because the scheme would be nonconservative, since the numerical flux F ·n at the common interface would have been computed in two different ways in Ω i and its neighbor and by using a q h that may already be non acceptable. Thus, the final P N P M solution in these cells, neighbors of troubled ones, is rearranged as follows: we keep the already computed values of the volume integral and of the surface integrals interacting with non-troubled cells, while the numerical flux across the troubled faces is substituted with the one computed through the limiter procedure.
Finally, note that for the subcell FV scheme we have a different CFL stability condition with CFL FV < 1 and h min the minimum cell size referred to Ω i . Condition (29) guides us in choosing the number of employed subcells N ω . In particular, our choice N ω = (2N + 1) d respects the stability condition (29) maintaining the original timestep size fixed for the current timestep (see (22)) but also taking into account the maximum possible number of subgrid elements allowed by that timestep size. We would like to stress again that our interest for P N P M schemes is motivated by the high resolution that they are able to provide and the reduced cost offered by the possibility of representing data with lower order polynomials of degree N and to achieve, however, an order of accuracy M , with M > N , using a very compact stencil and a simple reconstruction procedure.
The presented a posteriori subcell FV limiter is applied only where it is needed by detecting spurious oscillations a posteriori and it is based on strong stability preserving FV schemes developed precisely for dealing with discontinuous solutions. Since FV schemes are less accurate than P N P M schemes with N > 0, the limiter is applied on a finer subgrid than the original main grid in order to avoid a loss of useful information.

Adaptive Cartesian mesh refinement
The last ingredient that further increases the resolution of the proposed approach is the possibility of activating an Adaptive Mesh Refinement (AMR) technique based on a cell-by-cell refinement approach; indeed, the combined action of our subcell limiter and of AMR allows a sharp detection of all discontinuities. For details we refer to [15,14,73,49,43,120,54,53,24,116,93] and we recall here just the main features. Our algorithm basically consists in, starting from the main grid (2), introducing successive refinement levels, in regions of particular interest according to a prescribed refinement criterion. In particular, we have to fix the following parameters: the maximum level of refinement max , typically chosen equal to 2 or 3 in our tests; the refinement factor r, governing the number of subcells that are generated according to where ∆x is the size of the cell at refinement level number along the x-direction, and similarly for the other directions; the refinement criterion that we base on oscillations of second derivatives, see [81]. In practice, we have to compute where the summation k,l is taken over the number of space dimension of the problem in order to include the cross term derivatives, the parameter ε = 0.01 acts as a filter preventing refinement in regions of small ripples, and the function Φ = Φ(Q), that could be any suitable indicator function of the conserved variables Q, in our test is chosen to be simply Φ(Q) = ρ. Finally, the numerical solution at the subcell level during a refinement step is obtained by a standard L 2 projection, while a reconstruction operator is employed to recover the solution on the main grid starting from the subcell level. Moreover, in order to simplify the reconstruction procedure, the grid is treated as locally uniform for each cell independent of its grid level , because the neighbors cells at a coarser level − 1 can be virtually refined in order to allow for the reconstruction procedure on locally uniform meshes detailed in section 2.3. We also note that our AMR algorithm is endowed with a timeaccurate local time stepping (LTS) feature, see [49] for details.

Numerical results
In this Section we present a large set of numerical test cases in order to show the accuracy, robustness and efficiency of the presented P N P M family of schemes equipped with the a posteriori subcell finite volume limiter.
In order to cover a wide variety of physical phenomena we have applied our schemes to three sets of equations of relevance in fluid-dynamical applications, namely the Euler equations of compressible hydrodynamics (HD), the magnetohydrodynamics equations (MHD), and the special relativistic magnetohydrodynamics equations (RMHD).
In particular, for any set of equations we have selected both a smooth test case, to show the order of convergence of our schemes (up to order six), and some problems containing strong discontinuities going from logically onedimensional Riemann problems to classical challenging two-dimensional benchmarks, such as the Sedov explosion problem, the Double Mach Reflection problem, the MHD rotor problem, the RMHD blast wave, as well as the MHD & RMHD Orszag-Tang vortex problems. The presence of discontinuities allows to prove the robustness and resolution of our a posteriori subcell limiting strategy.
Moreover, the results obtained with the intermediate P N P M schemes (i.e. N = 0 and M > N ) are compared with the pure DG approach (i.e. N = M ) in order to show their gain in terms of computational efficiency, while maintaining a similar resolution. We also compare numerical results on AMR meshes against results obtained on fine uniform Cartesian meshes, demonstrating both the robustness of our schemes on adaptive meshes and the obtained savings in computational time.

Euler equations of gasdynamics
The first set of hyperbolic equations that we consider is given by the homogeneous Euler equations of compressible gasdynamics that can be cast in form (1) by choosing The vector of conserved variables Q involves the fluid density ρ, the momentum density vector ρv = (ρu, ρv) and the total energy density ρE. The fluid pressure p is related to the conserved quantities Q using the equation of state for an ideal gas where γ is the ratio of specific heats so that the speed of sound takes the form c = γp ρ .

Isentropic vortex
First of all, in order to verify the order of convergence of the proposed P N P M schemes, we consider a smooth isentropic vortex flow according to [68]. The computational domain is given by the square Ω = [0, 10] × [0, 10] with periodic boundary conditions set everywhere. For the initial conditions we consider a homogeneous background field Q 0 = (ρ, u, v, p) = (1, 1, 1, 1) traveling with a constant velocity v c = (1, 1) and we superimpose on this field some perturbations for density and pressure of the following form with the temperature fluctuation δT = − (γ − 1)ε 2 8γπ 2 e 1−r 2 and the vortex strength ε = 5. The velocity field is also affected by the following perturbations The initial condition is thus given by Q = Q 0 + δQ. The exact solution Q e at the final time t f can be simply computed as the time-shifted initial condition, In Table 2, we report the convergence rates from second up to sixth order of accuracy for the smooth vortex test problem run on a sequence of successively refined meshes up to the final time t f = 1.0. The optimal order of accuracy is achieved for the hybrid schemes P N P M with M > N and for the pure DG schemes with N = M .

The Sod shock tube problem
The Sod shock tube problem in 2D can be seen as a multidimensional extension of the classical Sod test case in 1D [108], which allows to verify the robustness and the resolution capacity of the employed numerical method on a rarefaction wave, a contact discontinuity and a shock at the same time, indeed the three waves are originated by the discontinuous initial condition.
Here, we consider as computational domain a square of dimension [−1, 1] × [−1, 1] covered with a uniform mesh of 50 × 10 control volumes, and the initial condition is composed of two different states, separated by a discontinuity at The final time is chosen to be t f = 0.4, so that the shock wave does not cross the external boundary of the domain, where wall boundary conditions are set. The algorithm for the calculation of the exact solution of this Riemann problem is given in [107].
We have run this problem with two fourth order methods, namely the P 2 P 3 and the P 3 P 3 schemes, and two sixth order methods, namely the P 3 P 5 and the P 5 P 5 schemes, equipped with the a posteriori subcell TVD FV limiter and employing the Rusanov flux both in the main P N P M scheme and at the limiter stage. The agreement of our numerical results with the exact solution is perfect and the hybrid schemes (i.e. M < N ) are faster than the pure DG schemes (N=M), see Figure 4.
Moreover, in Figure 5, one can see that the limiter activates exactly where the shock discontinuity is located also when the adaptive mesh refinement technique is employed.

The Lax shock tube problem
The Lax shock tube problem, introduced for the first time in [78], is another classical benchmark for high order methods for the solution of the Euler equations. The computational domain is the square [−1, 1] × [−1, 1] and the initial condition is composed of two different states, separated by a discontinuity located at x d = 0 In this case, we have covered our computational domain with a very coarse mesh of 20 × 10 elements activating the AMR procedure with max = 2 levels of refinement and r = 3. In Figure 6, we present the results obtained with two sixth order methods, namely the P 3 P 5 and the P 5 P 5 schemes, used together with the HLLEM [51,38] numerical flux. Both the numerical results perfectly agree with the reference solution and the hybrid scheme is 2.5 times faster Fig. 4: Sod shock tube problem at t f = 0.4. Left: we compare our numerical results obtained with two fourth order methods, namely the P 2 P 3 and the P 3 P 3 with the exact solution. Note that the two schemes show a similar resolution but the P 2 P 3 is twice faster than the P 3 P 3 ,Right: we compare our numerical results obtained with two sixth order methods, namely the P 3 P 5 and the P 5 P 5 with the exact solution. Note again that the two schemes show a similar resolution but the P 3 P 5 is 2.33 times faster than the P 5 P 5 . Fig. 5: Sod shock tube problem at t f = 0.4 solved with our P 3 P 5 sixth order scheme. We show in red the cells on which the limiter is activated and in blue the unlimited cells. The left panel is obtained with an initial grid of 50 × 10 elements and max = 2 levels of refinement with r = 3 in a total computational time of 1248s. The right image is obtained by using a fine uniform grid of 450 × 90 elements corresponding to the finest AMR grid level in a total computational time of 3670s. In both the cases the limiter perfectly activates only where the shock wave is located. than the pure DG scheme. Also the limiter activates exactly only at the shock location and, due to the subcell resolution, it does not affect the quality of the profile which is sharply captured even on a very coarse mesh.

The Shu-Osher shock tube problem
The Shu-Osher problem was first introduced in [100] and it allows to check the capability of our new scheme to deal simultaneously with physical oscillations and shock waves appearing at the same time during the simulation. It consists of a one-dimensional Mach 3 shock front interacting with a sinusoidal density disturbance that generates a combination of discontinuities and smooth structures, whose entropy fluctuations are amplified when passing through the shock.
We discretize our computational domain Ω = [−5, 5] × [0, 1] with an AMR grid with 64 × 4 elements on the coarsest level and a maximum refinement level max = 2 with r = 3. The initial conditions are given by and the simulations run up to the final time t f = 1.8. For this test case, we employ the HLL Riemann solver, and the TVD a posteriori subcell finite volume limiter. We then compare the results obtained with two pure DG schemes P 3 P 3 and P 5 P 5 , and the hybrid schemes P 2 P 3 and P 3 P 5 which have, respectively, the same order of accuracy. As shown in Figure 7, the methods are accurate, and robust thanks to the employed limiter strategy, and the intermediate schemes P N P M with M > N are computational more efficient than pure DG schemes. Fig. 6: Lax shock tube problem at t f = 0.14, obtained with our sixth order schemes, namely the P 3 P 5 and P 5 P 5 schemes. Left: we draw the density profile on the z axis and we colour in red the cells on which the limiter is activated and in blue the unlimited cells. Right: we compare our numerical results with the exact solution. Fig. 7: Shu-Osher shock tube problem at t f = 1.8 solved with our fourth order scheme, namely the P 2 P 3 and P 3 P 3 schemes (first column), and our sixth order schemes, namely the P 3 P 5 and P 5 P 5 schemes (second column). In the first two rows we plot the value of the density on the z axis and we depict in red the cell where the limiter is activated. In the third row our results are compared with a reference solution obtained with a WENO FV scheme on a very fine mesh. Moreover the computational time required by the P N P M schemes with N < M (respectively 394s and 2007s) are significantly shorter than those required by the pure DG schemes (N = M ) (respectively 880s and 5050s); nevertheless all the results show an excellent agreement with the reference solution. Fig. 8: Sedov test problem. Comparison of numerical density obtained with a selection of our high order methods (left, middle) and density contours obtained with the hybrid P 2 P 3 third order scheme on an AMR grid. Table 3: Sedov problem. We report the total CPU time in seconds and the values of the density peak (that should be equal to 6) obtained with a selection of methods going from order 4 to 8. The results are ordered with respect to the density peak value. One can notice that the hybrid schemes (as the P 2 P 3 and the P 3 P 5 ) have accurate results at a lower computational cost.

Sedov problem
This test problem is widespread in literature [71,83,57] and it describes the evolution of a blast wave that is generated at the origin O = (x, y) = (0, 0) of the computational domain Ω(0) = [0, 1.2] × [0, 1.2]. An exact solution based on self-similarity arguments is available from [97] and the fluid is assumed to be an ideal gas with γ = 1.4, which is initially at rest and assigned with a uniform density ρ 0 = 1. The initial pressure is p 0 = 10 −6 everywhere except in the cell V or containing the origin O where it is given by being E tot = 0.244816 the the total energy concentrated at x = 0 and |V or | the total volume of V or . Fig. 9: Limited cells (red) and unlimited ones (blue) for the Sedov problem solved with a selection of high order methods, i.e. P 1 P 5 , P 3 P 5 , P 5 P 5 , P 1 P 4 , P 2 P 3 , P 7 P 7 . A part some spurious oscillations of the P N P M schemes with N = 1, the limiter activates exactly at the shock location.
We solve this numerical test with a selection of high order methods going from fourth to eighth order of accuracy on a mesh of 50 × 50 elements with and without AMR. When the adaptive mesh refinement is activated, we take max = 2 and r = 3. For all these test cases we employ the Rusanov flux, CF L = 0.9, and the WENO a posteriori subcell finite volume limiter.
We show the results on the activation of the limiter in Figure 9, and the obtained density profiles in Figure 8. Finally, we compare the performances of a selection of methods in Table 3, in particular we highlight the needed computational times versus their capability of capturing the high density peak. We can remark that the hybrid schemes (as the P 2 P 3 and the P 3 P 5 ) have accurate and robust results at a lower computational cost; also the combination with adaptive mesh refinement helps in increasing the accuracy keeping the computational cost low.

Double Mach Reflection
The double Mach reflection problem was first studied by Woodward and Colella in [117] from which we take the setup also for our test.
We consider a computational domain Ω = [0, 4]×[0, 1] covered with a coarse mesh of 72 × 24 elements where we activate the adaptive mesh refinement with max = 2 and r = 3, and we compare the behavior of two fifth order Fig. 10: Double Mach Reflection density contours obtained with two fifth order schemes, namely the P 2 P 4 and the P 4 P 4 schemes. We plot 30 equally spaced contour lines from 1.5 to 22.9705 as suggested in [98]. schemes, namely the hybrid P 2 P 4 scheme and the P 4 P 4 pure DG scheme. For all the simulations, we employ the Rusanov flux and the second order TVD a posteriori subcell finite volume limiter.
At the beginning of the computation a shock wave moving at Mach number 10 is positioned at (x, y) = (1/6, 0) with an angle of 60 • with respect to the x-axis and the initial pre-shock conditions (on the left of the shock) are given by a constant density equal to 1.4 and a constant pressure p = 1. At the bottom boundary we employ reflective boundary conditions for x > 1/6 where we suppose the presence of a wall, and the exact post-shock conditions for 0 ≤ x ≤ 1/6 to mimic an angled wedge. At the top boundary, the flow variables are set to describe the exact motion of the Mach 10 shock. Finally at the left and right boundaries we set inflow and outflow boundaries.
The obtained numerical results are shown in Figure 10 for the entire domain; we also plot a zoom in Figure 11 where, one can notice the roll up of the Mach stem due to Kelvin-Helmholtz instabilities.

Ideal MHD equations
Next, we consider the equations of ideal classical magnetohydrodynamics (MHD) which, with respect to the previous set of equations, take also into account the evolution of the magnetic field B. The vector of the conserved variables Q and the flux tensor F of the general form (1) are given by Here, B = (B x , B y , B z ) represents the magnetic field and p t = p + 1 8π B 2 is the total pressure. The hydrodynamic pressure is given by the equation of state used to close the system, thus System (39) requires an additional constraint on the divergence of the magnetic field to be satisfied, that is ∇ · B = 0.
Here, (39) includes one additional scalar PDE for the evolution of the variable ψ, which is needed to transport divergence errors outside the computational domain with an artificial divergence cleaning speed c h , see [88,31]. A similar approach is adopted in [54,20,16]. A more recent and more sophisticated methodology to fulfill this condition exactly at the discrete level also in the context of high order ADER WENO finite volume schemes on unstructured simplex meshes can be found in [7].   O6 P 0 P 5 P 1 P 5 P 2 P 5 P 3 P 5 P 4 P 5 P 5 P 5 7.1e-01 0. First, for the numerical convergence studies, we solve the vortex test problem proposed by Balsara in [6]. The computational domain is given by the box Ω = [0, 10] × [0, 10] with wall boundary conditions imposed everywhere. The initial condition can be written in terms of the vector of primitive variables V = (ρ, u, v, w, p, B x , B y , B z , Ψ ) T as V(x, 0) = (1, δu, δv, 0, 1 + δp, δB x , δB y , 0, 0) T ,

MHD vortex
with δv = (δu, δv, 0) T , δB = (δB x , δB y , 0) T and where e z = (0, 0, 1), r = (x − 5, y − 5, 0) and r = r = (x − 5) 2 + (y − 5) 2 . The divergence cleaning speed is chosen as c h = 2. The other parameters are q = 1 2 , κ = 1 and µ = √ 4π, according to [6]. In Table 4, we report the convergence rates from second up to sixth order of accuracy for the MHD vortex test problem run on a sequence of successively refined meshes up to the final time t f = 1.0. The optimal order of accuracy is achieved both in space and time both for the hybrid schemes P N P M with M > N and for the pure DG schemes with N = M .

MHD rotor problem
The MHD rotor problem is a classical benchmark for MHD that was first proposed by Balsara and Spicer in [9]. It consists of a rapidly rotating fluid of high density embedded in a fluid at rest with low density. Both fluids are subject to an initially constant magnetic field.
The rotor produces torsional Alfvén waves that are launched into the outer fluid at rest, resulting in a decrease of angular momentum of the spinning rotor. We consider as computational domain the square Ω = [−0.5, 0.5] × [−0.5, 0.5] and as initial condition we take the density inside a circle of radius r ≤ 0.1 equal to ρ = 10, while the density of the ambient fluid at rest is set to ρ = 1. The rotor has an angular velocity of ω = 10. The pressure is p = 1 and the magnetic field vector is set to B = (2.5, 0, 0) T in the entire domain. As proposed by Balsara and Spicer we apply a linear taper to the velocity and to the density in the range from 0.1 ≤ r ≤ 0.105 so that density and velocity match those of the ambient fluid at rest at a radius of r = 0.105. The speed for the hyperbolic divergence cleaning is set to c h = 8 and γ = 1.4 is used. Wall boundary conditions are applied everywhere.
We run this problem on a coarse mesh made of 45 × 45 elements activating the AMR procedure with max = 2 levels of refinement and r = 3, and for comparison we also employ a finer uniform mesh of 405 × 405 elements corresponding to the finest AMR grid level. In particular, we have employed the P 2 P 4 fifth order scheme and the P 4 P 5 sixth order scheme with the Rusanov numerical flux and our a posteriori subcell WENO FV limiter. In all the cases, we can observe a good agreement between the obtained numerical results and those available in the literature, see

MHD Orszag-Tang vortex
We consider now the the vortex system of Orszag and Tang [89,30,91] for the ideal MHD equations. We choose as computational domain the square Ω = [0, 2π] × [0, 2π] with periodic boundary conditions set everywhere; we cover it with a uniform grid of 128 × 128 elements. The initial condition written in terms of primitive variables are the following (ρ, u, v, w, p, B x , B y , B z ) = γ 2 , − sin(y), sin(x), 0, γ, − √ 4π sin(y), √ 4π sin(2x), 0 (44) with γ = 5/3. The divergence cleaning speed is set to c h = 2 and the final time of the simulation is taken to be t f = 3 as in [39]. We solve this test by employing three different fifth order schemes, namely the hybrid P 2 P 4 and P 3 P 4 schemes and the pure DG P 4 P 4 scheme, with the Rusanov numerical flux and equipped with our a posteriori subcell TVD finite volume limiter. The obtained numerical results and the cells on which the limiter is activated are presented in Figure 14. One can notice that the three methods produce similar results with a good qualitative agreement compared to the solutions provided in [39,7,120]; moreover, the hybrid schemes are computationally more efficient than the pure DG scheme.  We compare the results obtained with three fifth order schemes, namely the P 2 P 4 , P 3 P 4 , P 4 P 4 schemes. In the first column we depict the density contours, in the second column we depict the pressure contours and in the third column we depict the troubled cells in red and the unlimited cells in blue. One can notice that the three methods lead to similar results but with respect to P 4 P 4 the P 2 P 4 is 3.39 times faster and the P 3 P 4 is 1.60 time faster.

Special relativistic MHD equations
The system of equations of special relativistic magnetohydrodynamics (RMHD) is supposed to provide a sufficiently accurate description of the dynamics of those astrophysical plasma that move close to the speed of light and which are subject to electromagnetic forces that dominate over the gravitational forces. For example this is the case of high energy astrophysical phenomena like extragalactic jets [13], gamma-ray bursts [77] or magnetospheres of neutron stars [86].
For a more detailed description of this model and a review of the numerical methods used in its approximation we refer to [119] and the reference therein.
Here we briefly recall only the main terms appearing in the equations, which indeed can be written under the general hyperbolic form (1) by choosing where we have employed the classical tensor index notation based on the Einstein summation convention, which implies summation over two equal indices.
The conserved variables (D, S j , U, B j ) are related to the rest-mass density ρ, to the thermal pressure p, to the fluid velocity v i and to the magnetic field B i by D = ρW, where ε ijk is the spatial Levi-Civita tensor and δ ij is the Kronecker symbol. As usual in ideal MHD, the electric field is given by E = −v × B. The spatial tensor W i j in (45), representing the momentum flux density, is where δ ij is the Kronecker delta. The above equations include the divergence free condition ∇· B = 0 for the magnetic field, which, although is guaranteed by the Maxwell equations at a continuous level, is not automatically satisfied from a numerical point of view. Different strategies can be adopted in order to solve this problem (see [111] for a review). Here, as for the MHD case of Section 3.2, we have adopted the so called divergence-cleaning approach presented in [88,31], which considers an augmented system with an additional equation for a scalar field Φ, in order to propagate away the deviations from ∇ · B = 0 while the fluxes for the evolution of the magnetic field are also modified, namely

Alfvén wave
As for the previous set of equations we first of all check the convergence of our new numerical scheme. In the case of RMHD we can consider the propagation of a circularly polarized Alfven wave, for which an analytic solution can be computed, see [75,32].
As initial condition we impose the following profile for the magnetic field and the velocity field where B 0 is the uniform magnetic field along x, ρ = p = B 0 = η = 1, k is the wave number, while v A is the Alfven speed at which the wave propagates and γ = 5/3. For the computational domain, we consider the 2D square Ω = [0, 2π] × [0, 2π] with periodic boundary conditions set everywhere, and we run our simulation up to the final time t f = L/v A = 2π/v A corresponding to one period.
In Table 5, we report the L 2 norm of the errors between our numerical results and the analytical solution for the variable ρ. The convergence rates from second up to sixth order of accuracy are confirmed both for the hybrid schemes P N P M with M > N and for the pure DG schemes with N = M .

Riemann problems
Now, in order to check the robustness and accuracy of our a posteriori subcell FV limiter for the general class of the P N P M schemes, we solve two classical Riemann problems of RHD (i.e. RMHD with B = 0) for which also an exact solution is available. We consider the computational domain Ω = [−0.5, 0.5]×[0, 1] and as initial condition we impose the discontinuous values given in Table 6. We solve these two test cases with a fifth order P 3 P 5 scheme and the HLLEM numerical flux, over a mesh of 20 × 10 elements with max = 2 levels of refinement and r = 3.
The obtained numerical results, see Figure 15, show once again that our limiter procedure preserves the resolution of the underlying P N P M scheme even on a coarse mesh.

Cylindrical blast wave
We now take into account a truly two dimensional test in RMHD, i.e. the cylindrical expansion of a blast wave in a plasma with an initially uniform  Table 6 for the initial conditions. In the Figure, we compare our numerical results for the density ρ (squares) with the exact solution (continuous line). Fig. 16: RMHD blast wave at time t f = 4.0. We show the results obtained with two sixth order schemes, namely the P 3 P 5 and P 5 P 5 schemes. In the left column we plot the density contours and in the right column we depict the troubled cells in red and the unlimited cells in blue. The resolution and the number of limited cells are quite similar with the two approaches but the hybrid P 3 P 5 scheme is 2.79 times faster than the pure DG P 5 P 5 scheme. magnetic field. This is a severe test proposed in [76], and subsequently also solved in [80,32,48,119].
For the initial condition we set the rest-mass density and the pressure equal to ρ = 0.01 and p = 1 within a cylinder of radius r = 1.0, and ρ = 10 −4 and p = 5 × 10 −4 outside. Like in [76] and in [32], the inner and outer values are joined through a smooth ramp function between r = 0.8 and r = 1, to avoid a sharp discontinuity in the initial conditions. The plasma is initially at rest and subject to a constant magnetic field along the x-direction, i.e. B x = 0.1, B y = 0, B z = 0.
We have solved this problem on the computational domain Ω = [−6, 6] × [−6, 6], with a uniform mesh of 160 × 160 elements. We have used the Rusanov numerical flux and two sixth order schemes, namely the P 3 P 5 and the P 5 P 5 schemes, equipped with the robust a posteriori subcell second-order TVD FV limiter. The obtained numerical results, which agree with those available in the literature, are reported in Figure 16.

RMHD Orszag-Tang vortex
Finally, we have chosen the relativistic version of the well known Orszag-Tang vortex problem, proposed by [89], and adapted to the relativistic case in [48,119]. The computational domain is Ω = [0, 2π] × [0, 2π] and the initial Fig. 17: RMHD Orszag Tang at time t f = 5. We present the density contours (left column) and the limited cells (in red in the right column) obtained with two sixth order schemes, namely the P 3 P 5 and the P 5 P 5 schemes, on an adaptive mesh with 45 × 45 control volume on the coarsest level, max = 2 and r = 3.
with γ = 4/3. To solve this system we employ the sixth order hybrid scheme P 3 P 5 over a level zero mesh of 45 × 45 elements, activating the AMR feature with max = 2 levels of refinement and r = 3. The obtained numerical results are reported in Figure 17: once again we can notice that the two schemes have a similar resolution but the hybrid scheme is 2.62 times faster than a pure DG scheme, of the same order. Furthermore, we can observe that the proposed a posteriori subcell limiter procedure is robust and maintains the high resolution of the underlying P N P M scheme even on coarse meshes.

Conclusion
In this paper we have proposed a new simple, robust, accurate and computationally efficient limiting strategy for the general family of ADER P N P M schemes, allowing, for the first time in literature, the use of hybrid reconstructed methods (N > 0, M > N ) in the modeling of discontinuous phenomena. The key ideas behind our limiter are: i) its local activation only where the linear schemes introduces oscillations through an a posteriori detector, ii) its robustness due to the use of strong stability preserving TVD or WENO FV schemes as limiter, iii) its resolution due to the use of the limiter on 2N + 1 subcells. Thus, we have been able to apply this new approach to many different systems of hyperbolic conservation laws, providing highly accurate numerical results in all cases. Moreover, we had the possibility to compare the performance of the class of intermediate P N P M schemes with M > N > 0 with pure DG schemes (M = N ). We have observed that in most cases the intermediate P N P M schemes offer a similar resolution compared to pure DG methods, but at a reduced computational cost.
Future work will consider the extension of P N P M scheme with N > 0, M > N to unstructured moving meshes [18,60], in particular for regenerating Voronoi tessellations [57,56], and to the three-dimensional case. Finally, due to their low memory consumption and their gain in computational efficiency compared to DG schemes, they will be also considered for astrophysical applications [58,42,40] and the unified first order hyperbolic model of continuum mechanics proposed in [90,46,47,25], where a large number of conserved variables has to be discretized. Due to their accuracy and compact stencil in the future we also plan to use P N P M schemes with a posteriori subcell finite volume limiter in the context of hyperbolic reformulations of nonlinear dispersive systems and wave propagation problems, see e.g. [33,55,52,11,12].