A Posteriori Subcell Limiting of the Discontinuous Galerkin Finite Element Method for Hyperbolic Conservation Laws

A Posteriori Subcell Limiting of the Discontinuous Galerkin Finite Element Method for Hyperbolic Conservation Laws

Michael Dumbser michael.dumbser@unitn.it Olindo Zanotti olindo.zanotti@unitn.it Raphaël Loubère raphael.loubere@math.univ-toulouse.fr Steven Diot diot@lanl.gov Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano, 77 - 38123 Trento, Italy. CNRS and Institut de Mathématiques de Toulouse (IMT) Université Paul-Sabatier, Toulouse, France Fluid Dynamics and Solid Mechanics (T-3), Los Alamos National Laboratory, NM 87545 U.S.A.
Abstract

The purpose of this work is to propose a novel a posteriori finite volume subcell limiter technique for the Discontinuous Galerkin finite element method for nonlinear systems of hyperbolic conservation laws in multiple space dimensions that works well for arbitrary high order of accuracy in space and time and that does not destroy the natural subcell resolution properties of the DG method. High order time discretization is achieved via a one-step ADER approach that uses a local space-time discontinuous Galerkin predictor method to evolve the data locally in time within each cell.

Our new limiting strategy is based on the so-called MOOD paradigm, which a posteriori verifies the validity of a discrete candidate solution against physical and numerical detection criteria after each time step. Here, we employ a relaxed discrete maximum principle in the sense of piecewise polynomials and the positivity of the numerical solution as detection criteria. Within the DG scheme on the main grid, the discrete solution is represented by piecewise polynomials of degree . For those troubled cells that need limiting, our new limiter approach recomputes the discrete solution by scattering the DG polynomials at the previous time step onto a set of finite volume subcells per space dimension. A robust but accurate ADER-WENO finite volume scheme then updates the subcell averages of the conservative variables within the detected troubled cells. The recomputed subcell averages are subsequently gathered back into high order cell-centered DG polynomials on the main grid via a subgrid reconstruction operator. The choice of subcells is optimal since it allows to match the maximum admissible time step of the finite volume scheme on the subgrid with the maximum admissible time step of the DG scheme on the main grid, minimizing at the same time the local truncation error of the subcell finite volume scheme. It furthermore provides an excellent subcell resolution of discontinuities.

Our new approach is therefore radically different from classical DG limiters, where the limiter is using TVB or (H)WENO reconstruction based on the discrete solution of the DG scheme on the main grid at the new time level. In our case, the discrete solution is recomputed within the troubled cells using a different and more robust numerical scheme on a subgrid level.

We illustrate the performance of the new a posteriori subcell ADER-WENO finite volume limiter approach for very high order DG methods via the simulation of numerous test cases run on Cartesian grids in two and three space dimensions, using DG schemes of up to tenth order of accuracy in space and time (). The method is also able to run on massively parallel large scale supercomputing infrastructure, which is shown via one 3D test problem that uses 10 billion space-time degrees of freedom per time step.

keywords:
Arbitrary high-order Discontinuous Galerkin schemes, a posteriori subcell finite volume limiter, MOOD paradigm, ADER-DG, ADER-WENO, high performance computing (HPC), hyperbolic conservation laws
journal: Journal of Computational Physics

1 Introduction

The discontinuous Galerkin (DG) finite element method has been originally proposed by Reed and Hill reed . Later, a solid theoretical framework has been established by Cockburn and Shu in a well-known series of papers cbs0 ; cbs1 ; cbs2 ; cbs3 ; cbs4 for the application of discontinuous Galerkin schemes to nonlinear hyperbolic systems of conservation laws. A very important property of DG schemes is that they satisfy a local cell entropy inequality for any polynomial degree used for the approximation of the discrete solution. As a consequence, this guarantees nonlinear stability in norm for arbitrary high order of accuracy, see the proof by Jiang and Shu jiangshu for the scalar case and its subsequent extensions to systems BarthCharrier ; HouLiu . This means that the DG scheme is by nature very robust and clearly appropriate for the solution of nonlinear hyperbolic conservation laws. However, being a linear scheme in the sense of Godunov godunov , even the DG method needs some sort of nonlinear limiting to avoid the Gibbs phenomenon at shock waves or other discontinuities. There is a vast literature on the topic of limiters for DG schemes and a non-exhaustive review on this topic will be presented later in section 3 of this paper. The key idea of many DG limiters is the following: first an unlimited solution is computed with the DG scheme, then an indicator detects so-called troubled cells, i.e. those zones of the domain which may need limiting, see Qiu_2005 for a detailed comparison. For troubled cells, the degrees of freedom of the discrete solution are then modified by some sort of nonlinear reconstruction technique, based on the discrete solution in the troubled cell and its neighbors.

Concerning the time discretization, mostly explicit TVD Runge-Kutta schemes are used, which lead to the so-called Runge-Kutta DG schemes. A review on DG schemes can be found in CBS-book ; CBS-convection-dominated . However, explicit DG schemes suffer from a very severe time step restriction where the maximum admissible Courant number typically scales as approximately , if denotes the polynomial degree of the approximation of the DG scheme. Alternative high order accurate explicit time discretizations for DG schemes have been explored in QiuDumbserShu and taube_jsc ; dumbser_jsc , which, however, have led to an even more restrictive CFL condition. While the DG method is mostly used only for spatial discretization, it has been introduced as a uniform discretization of space and time in the global space-time DG scheme of Van der Vegt et al. spacetimedg1 ; spacetimedg2 ; KlaijVanDerVegt , which leads to an implicit method of theoretically arbitrary high order of accuracy in space and time and which is unconditionally stable. The strategy presented in spacetimedg1 ; spacetimedg2 , however, requires the solution of a global nonlinear algebraic system at each time step. In order to reduce the complexity of globally implicit space-time DG schemes, in DumbserEnauxToro ; HidalgoDumbser ; Dumbser2008 a local space-time DG approach has been suggested, which leads only to an element-local implicit method. However, also in this case the final DG scheme is explicit and thus has to satisfy the typical stability condition of explicit DG schemes.

In the finite volume context, recently a new concept has been proposed, namely the Multi-dimensional Optimal Order Detection (MOOD) approach, which is an a posteriori approach to the problem of limiting. The key idea of this paradigm is to run a spatially unlimited high-order finite volume scheme in order to produce a so-called candidate solution. Then the validity of this candidate solution is tested against a set of predefined admissibility criteria. Some cells are marked as ’acceptable’ and are therefore valid. Some others may be locally marked as ’problematic’ or ’troubled’, if they do not pass the detection process. These cells and their neighbors are consequently locally recomputed using polynomial reconstructions of a lower degree. Thus, after decrementing the polynomial degree and locally recomputing the solution, a new candidate solution is obtained. That solution is again tested for validity and the decrementing procedure re-applies, if necessary.

Such order decrementing can occur several times within one time step for the same cell, but it will always halt after a finite number of steps: either the cell is valid for a polynomial degree greater than , or the degree zero is reached. In the worst case, a cell is updated with a robust and stable first order accurate Godunov-type finite volume scheme, which is supposed to produce always valid (monotone and positivity-preserving) solutions under CFL condition. This a posteriori check and order decrementing loop is called the ’MOOD loop’. We refer the reader to CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 for more details.

The link between the MOOD concept developed in the finite volume framework and the typical strategy adopted in a classical DG limiter seems obvious. In both approaches first a candidate solution is computed using an unlimited scheme. Then, troubled cells are detected based on some criteria and the discrete solution is corrected. However, there are important differences: the typical DG limiters postprocess the candidate solution using a nonlinear reconstruction technique. They furthermore use only one time level (the current one) for the detection and postprocessing step. In contrast, the MOOD approach first of all uses two time levels (the old one and the current one) for the detection of troubled cells. Second, the MOOD approach recomputes the solution using a different numerical scheme, which is supposed to be more robust at shock waves. The fact of recomputing the solution and of looking at two different time levels for the detection of troubled cells makes it in the notation of the authors of CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 an a posteriori approach.

The key innovation of the present paper is now the use of the aforementioned a posteriori MOOD paradigm as a limiter for high order DG schemes. Since simple order decrementing as in the MOOD approach for finite volume schemes would obviously destroy the natural subcell resolution capability of DG, a more sophisticated strategy is needed here. We therefore suggest to recompute the solution of troubled cells on a finer subgrid inside each cell, using a more robust but still very accurate one-step ADER-WENO finite volume scheme titarevtoro ; Balsara2013934 ; AMR3DCL . The data can be scattered from the main grid to the subgrid and gathered back via appropriate subcell projection and subcell reconstruction operators. Such operators are in principle well-known from high order finite volume schemes and spectral finite volume schemes spectralfv2d ; spectralfv3d ; spectralfv.dg . The choice of the subgrid size is very important. In this paper we suggest to choose the subgrid size so that the local CFL number on the subgrid is as large as possible, hence the maximum admissible time step size of the finite volume scheme on the subgrid matches the maximum admissible time step of the DG scheme on the main grid. This leads to subcells per space dimension for a DG scheme using a piecewise polynomial approximation of degree . For alternative subcell methods used as limiter for the DG method see CasoniHuerta2 ; Sonntag . However, none of these uses our a posteriori detection concept, nor do they recompute the solution in troubled cells via a better than second order accurate subcell finite volume method. Furthermore, the method proposed in Sonntag does not use the optimal subgrid size that allows to get the maximum admissible CFL number on the subgrid.

The rest of this paper is organized as follows. The one-step discontinuous Galerkin scheme used in this paper (ADER-DG) is presented in section 2. In section 3 we give a thorough but still non-exhaustive overview of existing a priori limiters for the DG method, such as slope/moment reduction, artificial viscosity or WENO-like techniques. Then in section 4 we present all details of our new subcell-based a posteriori approach. Next, section 5 gathers all numerical results for a large set of different test cases in order to assess the validity and robustness of our new a posteriori subgrid limiter. Smooth and non-smooth test cases are simulated. Finally, conclusions and perspectives are drawn in section 6.

2 The ADER Discontinuous Galerkin Method

In this paper we consider nonlinear systems of hyperbolic conservation laws in multiple space dimensions of the form

(1)

with appropriate initial and boundary conditions

(2)

where is the state vector of conserved quantities, and is a non-linear flux tensor that depends on state . denotes the computational domain in space dimensions whereas is the space of physically admissible states, also called state space or phase-space.
We solve this system of equations by applying the general family of methods introduced in Dumbser2008 , which provides high order of accuracy in both space and time. However, in this paper we only focus on the family of pure Discontinuous Galerkin (DG) schemes, i.e. the schemes in the context of Dumbser2008 . The numerical method is formulated as a one-step predictor corrector method GassnerDumbserMunz : in the predictor step (1) is solved within each element in the small (see eno ) by means of a locally implicit space-time discontinuous Galerkin scheme. The final time update of the discrete solution is explicit and is obtained by the one-step corrector. In the following we only summarize the main steps, while for more details the reader is referred to Dumbser2008 ; DumbserZanotti ; HidalgoDumbser ; GassnerDumbserMunz ; Balsara2013934 .

2.1 Space discretization and data representation

The computational domain is discretized by a Cartesian grid of conforming elements , where the index ranges from 1 to the total number of elements . The elements are chosen to be quadrilaterals in 2D and hexahedrons in 3D. The union of all elements represents the computational grid or the main grid of the domain,

(3)

We denote the cell volume by . At the beginning of each time-step, the numerical solution of equation (1) for the state vector is represented within each cell of the main grid by piecewise polynomials of maximum degree and is denoted by ,

(4)

where is referred to as the discrete “representation” of the solution. The space of piecewise polynomials up to degree is spanned by the basis functions . Throughout this paper we use the Lagrange interpolation polynomials passing through the tensor-product Gauss-Legendre quadrature points stroud as spatial basis functions, see also Kopriva ; GassnerKopriva ; KoprivaGassner .

2.2 Local space-time predictor

The representation polynomials are now evolved in time according to a local weak formulation of the governing PDE in space-time, see DumbserEnauxToro ; Dumbser2008 ; HidalgoDumbser ; DumbserZanotti ; GassnerDumbserMunz ; Balsara2013934 . The local space-time Galerkin method is only used for the construction of an element-local predictor solution of the PDE neglecting the influence of the neighbors. This predictor solution is further inserted into a corrector step described in the next section, which then provides the appropriate coupling between neighbor elements via a numerical flux function (Riemann solver).

Let us first transform the PDE (1) into a space-time reference coordinate system of the space-time reference element , with . The spatial reference elements are denoted by and are defined as , i.e. the unit square in two space dimensions and the unit cube in three space dimensions, respectively. Time is transformed according to . This yields

(5)

with the modified flux

(6)

To simplify notation, let us define the following two operators

(7)
(8)

which denote the scalar products of and over the space-time reference element and over the spatial reference element at time , respectively.

Now we multiply (5) with a space-time test function and subsequently integrate over the space-time reference control volume to obtain the weak formulation

(9)

The discrete solution of equation (9) in space and time, denoted by , as well as the discrete space-time representation of the flux tensor are assumed to have the following form

(10)
(11)

where we have used the Einstein summation convention over two repeated indices. If one employs a nodal basis as in Dumbser2008 , one simply has

(12)

Throughout this paper we use space-time basis functions that are the Lagrange interpolation polynomials passing through the tensor-product Gauss-Legendre quadrature points on the space-time reference element , see AMR3DCL . In that way, all the resulting matrices have a sparse block structure and the computations can be done efficiently in a dimension-by-dimension fashion. Integration by parts in time of the first term in (9) yields

(13)

Note that the piecewise high order polynomial representation is taken into account as initial condition of the Cauchy problem in the small in a weak sense by the term .
Eqn. (10) is then substituted into (13) and yields the following iterative scheme, see Dumbser2008 ; HidalgoDumbser ; DumbserZanotti for more details,

(14)

where is the iteration index. The iterative method (14) converges very efficiently to the unknown expansion coefficients of the local space-time predictor solution, see Dumbser2008 ; HidalgoDumbser ; DumbserZanotti . Once the are obtained the local space-time predictor is known inside each cell (10). The above iterative procedure has replaced the cumbersome Cauchy-Kovalewski procedure that has been initially employed in the original version of the ADER finite volume and ADER discontinuous Galerkin schemes schwartzkopff ; toro3 ; toro4 ; titarevtoro ; dumbser_jsc ; taube_jsc ; DumbserKaeser07 .

2.3 Fully discrete one-step ADER-DG scheme

The fully discrete one-step ADER-DG scheme dumbser_jsc ; taube_jsc ; QiuDumbserShu is obtained after multiplication of the governing PDE (1) by a test function , which is identical with the spatial basis functions, and subsequent integration over the space-time control volume . The flux divergence term is then integrated by parts and one obtains the weak formulation

(15)

where is the outward pointing unit normal vector on the surface of element . As usual in the DG finite element framework cbs0 ; cbs1 ; cbs2 ; cbs3 ; cbs4 the boundary flux term in (15) is then replaced by a numerical flux function (Riemann solver) in normal direction, , which is a function of the left and right boundary-extrapolated data, and , respectively. Inserting the local space-time predictor into (15) then yields the following arbitrary high order accurate one-step Discontinuous Galerkin (ADER-DG) scheme:

(16)

In general we use the simple Rusanov (local Lax Friedrichs) flux Rusanov:1961a , or the Osher-type flux recently proposed in OsherUniversal as numerical flux function at the element boundaries, although any other kind of Riemann solver could be also considered, see toro-book for an overview of state-of-the-art Riemann solvers.

2.4 Time step restriction

The DG method applied to convective problems enjoys a Courant-Friedrichs-Lewy (CFL) number that decreases with the approximation order , roughly it follows , CBS-book ; CBS-convection-dominated . This rather restrictive condition is caused by the growth of the spectrum of the spatial discretization operator of the semi-discrete scheme krivodonova2013analysis , hence the time step of an explicit DG scheme in multiple space dimensions has to satisfy

(17)

where and are a characteristic mesh size and the maximum signal velocity, respectively. This rather restrictive condition on the time step is to be put in perspective with the subcell resolution capability of any DG method, which allows the use of coarse or even very coarse grids. Recall that the representation polynomials are in principle of arbitrary degree . Consequently, within one cell multiple subcell features can be captured. We refer the reader to the next section for a detailed description about subcell resolution.

This closes the brief description of the family of one-step ADER Discontinuous Galerkin schemes that is further used in this paper. Note that this family of schemes is the unlimited version which, as such, can not prevent numerical oscillations and/or Gibbs phenomenon from occurring in the presence of steep gradients or shock waves. For sufficiently smooth solutions the unlimited ADER-DG scheme has shown an effective order of accuracy of , linked to the representation space . See Dumbser2008 for convergence tests in 2D and 3D. Nonetheless, adding some sort of limiting is of paramount importance to ensure the numerical stability of the method, despite the nonlinear stability properties that can be proven for the DG scheme in norm and the associated local cell entropy inequality jiangshu . It is the purpose of the next two sections to discuss the issue of limiters of DG schemes and to propose a new a posteriori subcell-based finite volume limiting strategy.

3 Nonlinear stability via a priori limiting

In this section we briefly review existing limiters for DG that have been designed upon paradigms such as slope limiters, artificial viscosity or essentially non oscillatory reconstruction procedures, such as ENO, WENO or HWENO. A discussion follows that will give birth to the paradigms and design principles of our a posteriori subcell-based MOOD limiting presented in the next section.

3.1 A brief and non-exhaustive review of a priori limiters for DG schemes

The main difficulty of solving nonlinear hyperbolic systems of conservation laws arises due to the fact that solutions of the system may become discontinuous, even if the initial conditions are smooth. This important discovery was due to the groundbreaking work of Bernhard Riemann riemann1 ; riemann2 . Despite their provable nonlinear stability jiangshu , even DG schemes can fail in the presence of strong shock waves or steep gradients and may generate strong oscillations that can ultimately lead to a failure of the computation. This is a consequence of the Godunov theorem godunov that also applies to the DG method. Therefore, some sort of non-linear limiting is needed and a vast literature on designing such techniques for DG methods does exist. Being exhaustive is almost impossible. As a consequence we only recall the underlying basic principles and some of the most important limiters used. Philosophically speaking, any limiter procedure for DG answers the following questions:

  • Q1: Where are the locations where limiting is needed? This is the so-called troubled zone indicator. We emphasize that the number of these locations may be very small, since for high order DG schemes many features can be captured within one cell length, see Fig. 1.

  • Q2: How do we achieve high order of accuracy along with a non-oscillatory property close to these locations? Ideally one should manipulate or replace the DG polynomials in such a way that additional numerical dissipation is supplemented close to the detected locations, but nowhere else, and preferably without destroying the subcell resolution of the DG method.

The design of such a priori limiters is a difficult task. However, amongst all DG limiters, three families seem to emerge: artificial viscosity based limiters such as Hartman_02 ; Persson_06 ; Barter_2010 ; Nguyen11 ; spacetimedg1 ; spacetimedg2 ; Feistauer4 ; Feistauer5 ; Feistauer6 ; Feistauer7 , “Slope” or moment limiters, for instance the total variation bounded (TVB) limiter of Cockburn and Shu cbs0 ; cbs1 ; cbs2 ; cbs3 , the moment-based limiters Biswas_94 ; Burbeau_2001 ; Yang_AIAA_2009 ; Yang_parameterfree_09 , or the hierachical slope limiters Kuzmin2010 ; Kuzmin2013 ; Kuzmin2014 , as well as WENO and HWENO based limiters, such as the ones developed in QiuShu1 ; QiuShu2 ; QiuShu3 ; Qiu_2005 ; balsara2007 ; Luo_2007 ; Zhu_2008 ; Zhong_2013 ; Zhu_3D_12 . Let us briefly describe and comment the design principles of these three families.

3.1.1 Artificial viscosity (AV)

The idea of using an artificial viscosity (AV) to stabilize shock waves dates back to von Neumann and Richtmyer vonneumann50 during the Manhattan project in the 1940’s at Los Alamos National Laboratory. This idea has been fruitful especially in the Lagrangian community since then and a lot of AV models emerged, bulk AV, edge AV, tensor AV, see pages 17-21 of Loubere_HDR_2013 and the reference herein for a review. Back to the origin, the illuminating idea of von Neumann and Richtmyer was to introduce a purely artificial dissipative mechanism of such a form and strength that the shock transition would be a smooth one, extending over a small number of cell lengths, and then to include this dissipation into the finite difference equations, see neu , page 312. A lot of minor and/or major improvements have been developed since then but any artificial viscosity technique revolves around the basic ideas: (i) define the region of the shock waves, (ii) add some dissipative mechanism over a small number of cells.
While it is nowadays rarely used in the context of shock capturing finite volume schemes, the artificial viscosity concept has become popular again in the context of DG schemes in order to capture shocks. In Hartman_02 the authors use the magnitude of the residual to determine the amount of artificial viscosity added to the shock region. In Persson_06 the authors have introduced a subcell shock-capturing method based on the manipulation of the density variable for determining the shock region and also the magnitude of artificial viscosity to be added. Piecewise constant artificial viscosity leading to oscillations has later been overcome in Barter_2010 where the AV model is based on the system of equations solved (one additional equation is solved for the AV). In the DG community researchers have improved the AV models adding artificial terms to the physical viscosity coefficients Cook:2004 or using an analytic function of the dilatation Nguyen11 111 This phenomenon also occurred in the Lagrangian community. For instance, the function of the dilatation in Nguyen11 is related to the so-called compression switch used in the Lagrangian community known from an unpublished work by Rosenbluth from Los Alamos in the 1950’s, where he suggested that the artificial viscosity should be zero when the fluid is undergoing an expansion. This ’trick’ is nowadays known as the ’artificial viscosity switch’ Csw98 ; Caramana-Burton-Shashkov-Whalen-98 ; Cam01 ..

3.1.2 WENO limiting procedure

Some other authors suggest the use of a (H)WENO limiting procedure for Runge-Kutta DG methods, for instance in QiuShu1 ; QiuShu2 ; QiuShu3 ; balsara2007 ; Zhong_2013 ; Zhu_13 . Usually these authors adopt the following framework:

  1. Identify invalid cells (called “troubled cells”), namely, the cells which might need limiting. In Qiu_2005 a very detailed comparison among different possible troubled-cell indicators is carried out. These indicators are a very important ingredient for triggering the subsequent WENO limiter procedure and are often based on minmod-type slope limiters, such as the modified TVB minmod limiter for the 2D scalar case cbs3 , the modified TVB minmod limiter in characteristic variables proposed in cbs2 and cbs4 for nonlinear systems in one and multiple space dimensions, respectively. Sometimes, other shock detection techniques such as the one suggested in Krivodonova2004 are used as troubled cell indicator. A review can be found in wang2011adaptive , chapter 6. A troubled zone indicator based on subcell information has been forwarded by Balsara et al. in balsara2007 .

  2. Replace the DG polynomials in these detected cells with reconstructed polynomials, that keep the original cell averages, preserve the same high order of accuracy, and are less oscillatory through a more or less classical WENO reconstruction procedure, see Zhong_2013 for the details. We underline that in the case of WENO limiters, the subcell resolution property of the DG method is lost to some extent, since the higher order moments are reconstructed from the cell averages defined on the coarse main grid.

Note that in the same spirit Hermite WENO schemes have also been designed to be used as limiters for DG schemes, see e.g. Qiu_2004 ; Luo_2007 ; Zhu_2009 ; balsara2007 , to avoid the non-compactness of the WENO limiters, but they are following a similar idea.
Usually improving the identification of the troubled cells is a key point, as the procedure may rely on user-defined parameters, which are often problem dependent even if in Zhong_2013 the authors have shown a relative insensibility of the method to these parameters. Last, the extension to 3D and more general meshes are important key points, see Zhu_3D_12 . More recently the troubled cell indicator has also been used for adaptive methods (mesh refinement), (order enrichment), or (mesh motion) by refining the troubled cells and/or coarsening the others Zhu_hadap_2013 .

3.1.3 ”Slope”/moment reduction

Slope limiting or moment reduction techniques may permit to control the jumps of the DG polynomials by constraining or nullifying the high-order components in designated cells BarthJespersen ; Burbeau_2001 ; cbs4 ; Yang_AIAA_2009 ; Yang_parameterfree_09 . A well-designed slope limiter must filter out non-physical oscillations without sacrificing the order of accuracy at smooth extrema. To do so some authors rely on monotonicity-preserving limiters frequently combined with ad hoc smoothness or oscillation detection procedures. Originally, slope limiters were developed to constrain piecewise linear polynomial reconstructions based on a discrete maximum principle, therefore constructing “slope” limiters for higher order polynomials are difficult to build into this paradigm Krivodonova2004 ; Michoski2011 . Many classical limiters have been already used, for instance the minmod-based TVB limiter cbs1 , the family of moment limiters Biswas_94 ; Burbeau_2001 , as well as monotonicity preserving limiters SureshHuynh ; Rider_2001 . Most of these limiters succeed in controlling spurious numerical oscillations, but they may also present a tendency to degrade accuracy when mistakenly used in certain smooth regions of the solution. In Kuzmin2010 ; Kuzmin2013 ; Kuzmin2014 the authors have developed a so-called hierarchical slope limiter. Briefly this limiter constrains the derivatives of polynomials, written in the Taylor polynomial form, in order to eliminate under/overshoots measured at the vertices of the cell, see Kuzmin2014 . This vertex-based hierarchical slope limiter has some common features with the moment limiter from Krivodonova2004 and seems to preserve smooth extrema without using any troubled cell indicator. Recently this hierarchical slope reduction problem has also been recast into an optimization problem Kuzmin2014 .

3.2 Discussion

Most of the aforementioned shock capturing procedures for DG schemes have a common feature: they rely on the fact that spurious numerical oscillations can be detected and corrected in the discrete solution by looking at only one time step and usually without using the PDE. The a priori character is evident for artificial viscosity based approaches, where the shock detector and the magnitude of the artificial viscosity are chosen based on the current solution . The TVB, (H)WENO and moment limiters can be to some extent considered as a posteriori limiter techniques, since they modify the higher order moments in troubled cells at the end of each time step (or Runge-Kutta stage) after having used an unlimited version of the DG scheme. Nonetheless, the troubled cell indicator as well as the limiter usually only consider the current discrete solution and higher order moments are replaced by some sort of nonlinear data reconstruction based on the degrees of freedom of the current solution, without using the PDE.

In a DG method many features can be captured within one characteristic cell length. Certain limiters may therefore dissipate a lot of these subscale features, and only subtle ones would maintain at the same time the high subcell accuracy and assure stability of the numerical solution at shock waves. According to numerical evidence provided in CasoniHuerta1 for the one-dimensional case, the artificial viscosity approach Persson_06 seems to be more appropriate for capturing subcell features compared to other DG limiters, especially for very high polynomial degrees . Ultimately, one may demand that a DG limiter acts only on the smallest length scale within one cell to avoid that excessive numerical dissipation impacts all features that are represented by . This smallest length scale is related to the size of the cells and to the polynomial degree .

Recent developments have been made concerning the construction of sub-cell limiters, which use either a finite volume method on subcells Sonntag or a smooth switch between a high order DG scheme and a first order finite volume subgrid method, CasoniHuerta1 ; CasoniHuerta2 , but none of these is based on the a posteriori concept proposed in this paper, nor do they rely on the use of higher order finite volume schemes on the subgrid.

Following the ideas of a novel a posteriori detection approach (MOOD), that has been introduced for the first time in the context of very high order accurate finite volume schemes in CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 , we propose to extend the MOOD concept to the context of very high order accurate DG finite element schemes in the following.

4 Nonlinear stability via a posteriori sub-cell limiting (SCL)

4.1 General design principles

The subcell resolution is also evident from the amount of data that are stored per cell in the DG context, namely from the number of the degrees of freedom, which is a function of the polynomial degree . For our tensor-product basis functions, the number of spatial degrees of freedom per cell is .

As already discussed before, limiting acts in two steps: First the troubled cell indicator detects which regions of the computational domain need limiting and, second, the limiting effectively adds some sort of numerical dissipation in these regions, either directly via artificial viscosity or via a nonlinear data reconstruction or slope limiting procedure. However, an appropriate DG limiter should ideally detect and correct problematic situations on a subscale level.

Starting with piecewise high order polynomial data it seems difficult to detect a priori problematic situations which will occur between time and time . In CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 in the context of high order finite volume (FV) schemes the authors have adopted a different strategy, called MOOD. This latter consists in testing a posteriori a so-called candidate solution that has been obtained using an unlimited high order scheme. If this candidate solution does not fulfill a set of properties called detection criteria, then the cells for which it has failed are recomputed using a more robust and more dissipative lower order scheme (decrementing of the order of accuracy of the scheme). For an invalid cell, the iterative MOOD loop ends either with a valid solution that has passed the detection criteria, or, with a solution that is updated with the lowest order scheme that is supposed to be monotone and positivity preserving under CFL condition. The detection criteria developed in CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 for different hyperbolic systems of conservation laws have proven to be sensitive enough to avoid excessive numerical diffusion while maintaining an effective high order of accuracy for smooth problems. More important, they are sufficient to dissipate numerical oscillations and ensure stability.

4.2 Subcell data representation, projection and reconstruction

If is the data representation of the DG scheme within cell at time and we consider a fine subgrid of , denoted by made of subcells called , , with , then we introduce an alternative data representation denoted by , which is defined by a set of piecewise constant subcell averages . These subcell averages are directly computed from by projection, which in this case simply means the computation of the integral average of since is piecewise constant on the subcells :

(18)

The above equation (18) is in the following also called projection operator , denoted by . Throughout this paper, the subcells are chosen to be equidistant Cartesian subcells AMR3DCL . To gather back the piecewise constant subcell data into a high order DG polynomial, we apply the following reconstruction operator : find and therefore so that

(19)

This is a classical reconstruction or recovery problem of a higher order polynomial from known cell averages, which typically arises within the finite volume context and also in spectral finite volume methods spectralfv2d ; spectralfv3d ; spectralfv.bnd . The approach presented in this paper explicitly admits , hence the resulting system may be overdetermined. This overdetermined system of reconstruction equations is then solved using a constrained least-squares approach, see kaeserjcp ; DumbserKaeser06b , where the constraint is the integral conservation of the cell average over the big cell , i.e.

(20)

The operator given by the solution of (19) and (20) is in the following denoted by . It is obvious that finding from (19) corresponds to the computation of the (pseudo-) inverse of the matrix associated with the projection operator (18). Hence, the two operators and satisfy the property

(21)

where is the identity operator. If the subgrid is large enough () then there is no loss of accuracy when applying the projection operator, since can always be recovered identically from using the recovery operator due to relation (21), thus is able to represent all the information contained in , see Fig. 1. If , like in Sonntag , then the number of subcells corresponds exactly to the number of degrees of freedom associated with the space .

Consequently, it is equivalent in terms of nominal accuracy to represent data either with a DG scheme of polynomial degree on the cell , or, representing the data by piecewise constant cell averages on the subgrid .

The use of subgrid information has been used very successfully also in the context of semi-implicit finite volume methods for free surface flows, see Casulli2009 ; CasulliStelling2011 , where it has led to a significant improvements in terms of numerical accuracy and computational efficiency.

Figure 1: Examples of DG polynomials on a cell (red) and associated projection onto subcell averages (blue). The information contained in can be recovered from via the subcell reconstruction operator for . Throughout this paper we use subcells.

4.3 Extension of MOOD to DG schemes

The a posteriori MOOD concept is now extended to the DG context as follows. First, a candidate solution is computed from by the unlimited DG scheme (16). Next, we apply the following detection criteria.

Physical admissibility detection (PAD)

The detection criteria must contain physics-based admissibility properties, and, as such can not be uncorrelated with the hyperbolic system of conservation laws which is solved. Hence, a candidate solution is said to be physically valid inside cell for this system if

(22)

where is the physical quantity that must satisfy the positivity constraint and which is a function of the state vector . For the Euler equations of compressible gas dynamics, for example, the mass density and the fluid pressure must be positive everywhere and for all times, hence and . There is a growing interest in designing high order accurate and positivity preserving finite volume and DG methods. For the most recent developments in the field of DG schemes see, for example, ShuPositivity1 ; ShuPositivity2 ; ShuPositivity3 .

Numerical admissibility detection (NAD)

In the past, the discrete maximum principle (DMP) was a very successful guideline for the construction of high resolution shock capturing schemes. In this paper, we therefore use the following relaxed version of a discrete maximum principle, which takes into account the data representation of the DG method under the form of piecewise polynomials. It is a natural extension of the DMP for cell averages used in CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 . The DMP is applied in an a posteriori manner as follows. A candidate solution is said to fulfill the numerical admissibility detection criterion in cell if the following relation is fulfilled componentwise for all conserved variables:

(23)

where is a set containing element together with its Voronoi neighbor cells that share a common node with . We see from (23) that the discrete maximum principle is now applied in the sense of polynomials. The polynomial that represents the candidate solution on element must remain between the minimum and the maximum of the polynomials that have represented the discrete solution at the old time step in the neighborhood of cell .

The quantity is used to relax the strict maximum principle in order to allow some very small overshoots and undershoots and to avoid problems with roundoff errors that would occur when applying (25) in a strict way, without using . Throughout this paper we set

(24)

with . We underline that the slight relaxation of the maximum principle has no influence on the positivity of the solution, because the positivity is detected separately under the PAD above. Since it is not very practical from a computational point of view to calculate the maximum and the minimum of the discrete solution within the neighborhood exactly, we use the following discrete version of (23) on the subscale level instead, which can be easily evaluated on the basis of subcell averages:

(25)

with and given by (27). A candidate solution is said to be valid inside cell , if it has passed both the physical and the numerical admissibility detection criteria, i.e. if (22) and (25) are satisfied. In this case, we set a cell-based indicator function to . If a cell does not fulfill the above a posteriori MOOD detection criteria, a cell is marked or flagged like in a typical troubled zones indicator by setting . However, we emphasize again that in our approach the detector uses information from two different time levels, and , while a classical troubled zone indicator would only look at the discrete solution at one time level. An alternative self-adjusting a priori strategy to trigger the subcell finite volume limiter presented in this paper could be the flattener algorithm proposed in PositivityWENO2 .

After the detection phase given by (22) and (25), the next operation of our subcell limiter consists in updating the discrete solution in invalid cells using a more robust scheme on the subgrid and based on the alternative data representation . For this task, one could in principle choose a simple and cheap second order TVD finite volume scheme, as used in Sonntag , but we prefer to use a higher order one-step ADER-WENO finite volume method to avoid the clipping of local extrema on the subgrid. This is particularly important if local extrema on the subgrid are associated with physical phenomena, such as sound waves. The WENO scheme on the subcells is able to resolve the subscale features without sacrificing neither the high resolution of smooth features, nor the robustness at shock waves and other discontinuities.

For invalid cells the discrete solution is then recomputed starting from the alternative data representation , given by piecewise constant subcell averages. The update is carried out using a third order ADER-WENO finite volume scheme on the Cartesian subgrid. For details see titarevtoro ; Balsara2013934 ; AMR3DCL , where all implementation details for reconstruction and one-step time update are given. Here, we just abbreviate the entire ADER-WENO scheme acting on the subcell averages by

(26)

If the cell has been flagged also in the previous time step as a troubled cell, the initial data for are directly available on the subgrid from the ADER-WENO finite volume scheme of the previous time step, otherwise, the initial data are given by the projection operator applied to the DG polynomials on the main grid. Appropriate boundary conditions are also needed for the WENO reconstruction on the subgrid of cell , like in the context of high order cell-by-cell AMR schemes Mulet1 ; AMR3DCL ; GrothAMR , hence the neighbor cell data are also scattered onto virtual subcells in the same way. To summarize, the initial and boundary data for the subcell ADER-WENO finite volume scheme are given by

(27)

Note that this operation is local and involves only the cell and its direct neighborhood , hence our subcell limiter (26) together with the necessary initial and boundary conditions (27) fits well in the general philosophy of DG schemes. Note also, that both, the ADER-DG scheme (16) as well as the ADER-WENO subcell limiter (26) are one-step schemes, hence the limiter is applied only once per time step, without the need of carrying out the same procedure in each Runge-Kutta substage of a classical RK-DG scheme again.

Finally, for troubled cells () we gather back the subgrid data representation produced by the subcell limiter using the reconstruction operator , which computes the final representation of the high order DG polynomial of degree on the main grid, i.e. . This concludes the description of our extension of the MOOD paradigm to DG schemes.

We stress that in our a posteriori subcell limiter approach the new discrete solution is recomputed by using a different and more robust numerical scheme,

(28)

while the TVB and (H)WENO limiters post-process the unlimited candidate solution by a nonlinear reconstruction operator, acting on selected degrees of freedom of the cell and its neighbors. Our limiter approach presented here is also very different from the one proposed by Sonntag and Munz Sonntag , who use an a priori switch from a DG formulation on the main grid to a subcell TVD finite volume scheme based on some a priori indicator function. It is also very different from the subcell approach forwarded in CasoniHuerta2 , which smoothly switches between a high order DG scheme and a first order finite volume subcell method based again on an a priori indicator function.

Figure 2: Top: classical DG algorithm with a priori limiter embedded in the solver. After the fluxes and the limiter have been computed no more action is taken and the solution can not be corrected if needed anymore — Bottom: sketch of our approach, where an unlimited candidate solution is provided by the DG solver. Subsequently, the scatter step projects the candidate solution onto the subgrid to get for applying the subgrid detection criteria. An a posteriori MOOD-type procedure is further applied based on the relaxed discrete maximum principle (25) and some physical admissibility criteria like positivity (22). Valid cells are left unmodified, whereas bad cells are recomputed using an ADER-WENO finite volume scheme on the subgrid. Last, for these troubled cells, a gathering step reconstructs the piecewise polynomials of degree on the main grid. The green color corresponds to operations made on the subgrid, black and blue colors correspond to cell based operations on the main grid.

4.4 Summary of the a posteriori subcell based MOOD limiting

The a posteriori subcell based MOOD limiter is a five-act play, see Fig. 2. We start with the cell centered piecewise polynomial representation of degree on the main grid in each cell and perform the following steps:

  1. High-order cell-based ADER-DG update. From the unlimited ADER-DG scheme (16) is used to compute a candidate solution . This DG scheme is meant to be the highest order accurate scheme that one wishes to use.

  2. Cell-to-subcell scattering (projection ). Project onto the subcells, . The subgrid of any cell must be dense enough to recover the original degrees of freedom of via the subcell reconstruction operator , hence .

  3. A posteriori MOOD detection procedure. Given a set of detection criteria like positivity (22) and DMP (25), determine on a subcell basis which cells of the main grid contain an invalid candidate solution . Mark the associated cell as invalid by setting . Otherwise, if all detection criteria are satisfied on a subcell basis, the cell does not require limiting and we set the flag . In this case the algorithm can exit with the candidate solution := without further treatment.

  4. Subcell-based ADER-WENO update. For each invalid cell project the initial and boundary data onto the subcells using (27). Then update the subgrid data using a higher order ADER-WENO finite volume scheme on the subgrid (26) to obtain . Note that the initial and boundary conditions for the ADER-WENO scheme are provided in each grid cell and its neighbors either by the projection operator (for unlimited cells), or by the ADER-WENO scheme applied in the previous time step (for limited cells).

  5. Subcell-to-cell gathering (reconstruction ). For any troubled cell with gather the new subgrid information into a cell-centered DG polynomial of degree on the main grid by applying the subcell reconstruction operator (19), .

Note that in this paper we do not consider the full successive order decrementing loop proposed in the original references on MOOD schemes, CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 , but we consider only two possible schemes: either the unlimited high order ADER-DG scheme on the coarse grid, or the more robust ADER-WENO finite volume scheme on the subgrid. This means that each grid zone is decremented at most once. Note also that, of course, the numerical fluxes in unlimited DG cells adjacent to limited cells must be recomputed in order to maintain conservation and consistency of the scheme. Within the DG finite element framework and the finite volume framework used here, however, it is no problem to deal with such hanging nodes at nonconforming grid interfaces properly, see AMR3DCL ; AMR3DNC for details in the context of high order finite volume schemes on space-time adaptive meshes with hanging nodes in space and time. For the DG method, the flux integral over the element boundary in Eqn. (16) is simply written as the sum of integrals over the sub-edges between the main grid and the subgrid. On both grids and for both methods (ADER-DG and ADER-WENO finite volume scheme), a space-time predictor needed for one-step numerical flux integration is available.

To summarize: handling non-conforming grids with different mesh size and different orders of accuracy in each zone (so-called -refinement) is possible in a very natural manner in the combined DG and finite volume framework used here, see also CBS-book . However, it means that if a cell is flagged for a posteriori limiting, those direct neighbors of the cell which are not marked for limiting also need to be recomputed to ensure conservation and consistency. Note, however, that only direct neighbors are affected and that the same treatment is also necessary in the original MOOD method, CDL1 ; CDL2 ; CDL3 ; ADER_MOOD_14 . Furthermore, the recomputation affects only those edges of unlimited cells that are adjacent to a limited element, hence the simplest way of practical implementation is a flux-correction approach that subtracts the original unlimited DG fluxes and adds the new subcell ADER-WENO fluxes across the respective subedges to all moments of the unlimited DG cell. Of course, this correction step of unlimited zones adjacent to troubled zones needs additional MPI communication within a parallel implementation of the scheme, but only direct neighbors are involved, hence the scheme is still fully local.

Concerning the additional memory requirements of our algorithm we would like to emphasize that the alternative data representation needs to be stored only in the troubled cells, while for the detection criterion (25) it is sufficient to store for each cell the maximum and minimum of each conserved variable found in the subcells of the neighborhood . If the limiter is acting only in a restricted number of zones of the computational domain, which is usually the case, then the memory overhead produced by our approach is very small. Last but not least, the same a posteriori detection approach is also applied to the projection of the initial condition, i.e. if the DMP and the PAD are not satisfied for with respect to the exact initial condition , we immediately activate the subcell representation in troubled zones at the initial time .

4.5 On the optimal choice of the subgrid size

At a first glance, the natural choice for the number of subgrid cells seems to be , in order to represent exactly the same amount of information within the subcell averages as originally contained in the high order DG polynomials on the main grid. This choice has been made for example in Sonntag . However, we are convinced that this choice is not the best one in terms of accuracy and local truncation error of the finite volume scheme on the subgrid. For a very detailed modified equation analysis (differential approximation) of high order ADER finite volume schemes, together with their dissipation and dispersion properties, see dumbser_diffapprox. In dumbser_diffapprox it was shown for the linear scalar advection equation in 1D that the error terms of ADER finite volume schemes of up to order 16 in space and time contain the factor , which means that the schemes are the more accurate the larger the CFL number. Looking at the typical (severe) time step restriction of explicit DG schemes given by (17), we note that by using subcells the CFL number for the finite volume method on the subgrid is only about half of the maximum admissible CFL number of an explicit Godunov-type finite volume scheme, since the finite volume method on the subgrid must satisfy the stability condition

(29)

Note that is the size of the subgrid cells. If we want to use the optimal time step on the subgrid with a Courant number close to the maximal one, the subgrid must satisfy , which clearly follows from (17) and (29). A coarser subgrid would lead to small local CFL numbers for the finite volume scheme on the subgrid and thus to more numerical dissipation and dispersion, in addition to the reduced sub-cell resolution due to the coarser mesh itself! In other words, the optimal value of yields lower errors not only due to a smaller mesh size, but also in terms of smaller constants in front of the error terms in the local truncation error analysis. Of course, also the choice is suboptimal, since in this case the subgrid finite volume scheme would limit the (already small) time step of the DG scheme on the main grid. At the end, the choice of is left to the user, but the aim of this section was to discuss the choice of the optimal value of in terms of accuracy and computational efficiency.

5 Numerical results

In this paper we focus on the the Euler equations of compressible gas dynamics

(30)

where denotes the mass density, is the velocity vector, is the fluid pressure, is the total energy density and denotes the identity matrix. With the notation we intend the dyadic product of the velocity vector with itself. To close the system the equation of state (EOS) of a perfect gas with adiabatic index is used:

(31)

The a posteriori MOOD-type subcell finite volume limiter has been implemented within an MPI parallel 3D code devoted to hyperbolic system of conservation laws on Cartesian grids, see AMR3DCL for the general framework. For most test cases presented in this section we have employed the ADER Discontinuous Galerkin method with piecewise polynomials of degree or , referred to as ADER-DG- and ADER-DG-, in the following. This DG method is then supplemented by the new a posteriori sub-cell limiter (SCL). The scheme acting at the subcell level is a third order ADER-WENO finite volume method titarevtoro ; AMR3DCL with reconstruction, denoted by WENO3. Hence the full scheme is called DG-+WENO3 SCL. For comparison purposes we will also employ low order DG schemes with in one test case.

In the Discontinuous Galerkin finite element framework, each variable inside a computational cell is represented by a polynomial of degree . In what follows, we will therefore always use the following visualization technique for 2D and 3D contour plots. We plot the numerical solution as pointvalues interpolated onto a subgrid with if the cell is unlimited () and we will plot the alternative representation of the solution on the subgrid used for the limiter with if the cell has been detected as a troubled cell (). For 1D cuts we usually take equidistant samples at subgrid level of the solution representations and , respectively. This is extremely important in order to verify that the subscale structure of the DG polynomials really represents a physically valid state within one large cell. Moreover, an acceptable DG method should maintain smooth solutions and, more importantly, it should produce genuinely discontinuous profiles at shock waves without spurious oscillations. In order to visualize which cells have been limited, in the rest of this section we will systematically represent in blue the unlimited cells () and in red () the limited cells, see for instance Figure 4 below.
A well known difficulty of a high order accurate method dealing with discontinuous solutions is to deposit entropy (and dissipation) on a length scale which is much smaller than the characteristics length of the coarse cell. Indeed, if the coarse cell is very large, spreading a shock wave over one or two of such cell(s) typically generates an excessively large numerical dissipation. Our a posteriori subcell limiter is supposed to act differently and in this section we will provide the numerical evidence of this property.
The methodology of validation is based on a sequence of classical test cases, namely:

  • Sod and Lax shock tube toro-book - these shock tubes are classical tests to assess the ability of a numerical method to deal with simple waves (rarefaction, contact discontinuity and shock wave).

  • Smooth vortex - this test is designed to observe the high-order of accuracy of a numerical method when the exact solution is smooth. Note that, in theory, in this case any limiter should not be activated at all if the mesh is fine enough. A detailed numerical convergence table for polynomial degrees ranging from to is provided.

  • Shu-Osher oscillatory shock tube shuosher2 - this test is designed to show the difficulty of capturing smooth small-scale features and shock waves.

  • Double Mach reflection woodwardcol84 - this classical test is meant to measure visually the ability of a numerical method to capture complex patterns created after the interaction of shock waves. We use this test to show the behavior of the method when the polynomial degree of the DG scheme is increased.

  • Forward facing step woodwardcol84 - this classical test simulates a supersonic flow over a forward facing step. The solution approaches a steady-state solution composed of multiple interacting shock waves and vortex like structures which are often dissipated with low order accurate schemes.

  • 2D Riemann problems kurganovtadmor - this classical suite of test problems is meant to assess the ability of a numerical method to solve genuinely two-dimensional Riemann problems emerging from four piecewise constant states joining at the origin.

  • Shock-vortex interaction Shock_vortex_95 ; caishu93 - this test is designed to observe the interaction of a planar shock wave with a cylindrical vortex and the subsequent complex structures of primary and secondary waves.

  • 3D explosion problem toro-book - this test is designed to show the ability of a scheme to deal with separate spherical waves in 3D and subsequently validates the approach in 3D. We use this test also to prove that our approach is well suited for being used within a massively parallel MPI framework running on 8000 CPU cores and capable of dealing with 10 billion space-time degrees of freedom per time step and conserved variable.

5.1 Sod and Lax shock tube

Here, we run the planar Sod shock tube problem and the classical Lax shock tube problem on a 2D structured mesh to assess the ability of the methods to capture one-dimensional simple waves. The initial conditions for density, velocity component and pressure are listed in Table 1. The other velocity component is initialized with .

Problem Left state Right state Final time
Sod 1.0 0.0 1.0 0.125 0.0 0.1 0.2
Lax 0.445 0.698 3.528 0.5 0.0 0.571 0.14
Table 1: Initial left and right states for the density , velocity and the pressure for the Sod and Lax shock tube problems. Final simulation times are also provided.

The ratio of specific heats is and for both problems the initial discontinuity is located in . The exact solution for these one-dimensional Riemann problems can be found in toro-book . The computational domain is paved with a very coarse structured mesh made of cells with and , see Figure 4. Dirichlet boundary conditions are imposed in -direction while periodic boundaries are applied in -direction.
In Figure 3 we present the results for DG- where density, velocity in -direction and pressure are displayed along the direction versus the exact solution. In this figure each DG- polynomial is represented by sample points. We observe a very good agreement with the exact solution, the shock being resolved just in one cell. In Figure 4 we also represent the unlimited polynomials (blue) and the limited ones (red), we can clearly see that the contact discontinuity in the Sod problem is resolved within one single cell in an almost S-type shape, thanks to the use of very high order polynomials of degree . On the other hand, although the cells embracing the shock region have been limited (they are colored in red), the subcell data are first of all non-oscillatory, second, the plateaus before and after the shock are well captured and third the shock is properly represented by a jump at an element interface. Back to the 1D panels we can see that this jump on subcell values is perfectly located at the exact shock wave location. The same comments also hold for the Lax problem. In conclusion, although the number of cells is very small, the ADER-DG- scheme with SCL can produce a very accurate solution and also the discontinuities are very sharp even on a very coarse mesh. The a posteriori limiter has acted properly to maintain the shock wave within one or two subcells.

Figure 3: Sod shock tube problem (left panels) at and Lax problem (right panels) at . In both cases a very coarse mesh of only cells on the main grid has been used. An ADER-DG- scheme supplemented with a posteriori ADER-WENO3 subcell limiter has been used — 1D cut on 200 equidistant sample points through the numerical solution (symbols) vs exact solution for density (top), velocity component (middle) and pressure (bottom).
Figure 4: Sod problem (top panel ) and Lax problem (bottom panel ) solved on a element mesh using ADER-DG- with a posteriori WENO3 subcell limiter — The density variable is displayed. Troubled cells are shown in red, while blue cells are updated with the unlimited ADER-DG- on the main grid.

5.2 Smooth vortex

The isentropic vortex problem was initially introduced for the two-dimensional compressible Euler equations in Shu1 to test the accuracy of numerical methods, since the exact solution is smooth and has a simple analytical expression. Let us consider the computational domain and an ambient flow characterized by , , , , , with a normalized ambient temperature computed with the perfect gas equation of state and .
A vortex is centered on the axis line at and supplemented to the ambient gas at the initial time with the following conditions , , , where

with and . The vortex strength is given by and the initial density follows the relation

(32)

Periodic boundary conditions are prescribed everywhere, so that at the final time the vortex is back to its original position. The numerical flux used here was the Osher-type flux presented in OsherUniversal . This problem has a smooth solution and thus should be simulated with effective high-order of accuracy if the limiter behaves properly. Four successively refined grids made of squares are employed to compute these errors. We compute the discrete , and error norms between the exact solution and the numerical solution for the density at the final time. The computation of these error norms is performed using a sufficiently high order accurate Gaussian quadrature rule. The errors and the rate of convergence are reported in Table 2 for the ADER-DG- supplemented with the a posteriori WENO3 subcell limiter and varying from to . From this table we can observe that the optimal order of convergence is essentially achieved by the scheme and that the proposed a posteriori subcell limiter does not destroy the accuracy of the high order DG scheme. Beyond it seems more difficult to get perfect orders due to the fact that the errors are very small and roundoff starts to play a role. However, even on these ultra coarse meshes the results are very accurate, for instance comparing the ADER-DG- on cells which is roughly as accurate as ADER-DG- on cells or ADER-DG- on elements.

2D isentropic vortex problem — ADER-DG- + WENO3 SCL
error error error order order order Theor.
DG- 25 9.33E-03 2.07E-03 2.02E-03 3
50 6.70E-04 1.58E-04 1.66E-04 3.80 3.71 3.60
75 1.67E-04 4.07E-05 4.45E-05 3.43 3.35 3.25
100 6.74E-05 1.64E-05 1.82E-05 3.15 3.15 3.10
DG- 25 5.77E-04 9.42E-05 7.84E-05 4
50 2.75E-05 4.52E-06 4.09E-06 4.39 4.38 4.26
75 4.36E-06 7.89E-07 7.55E-07 4.55 4.30 4.17
100 1.21E-06 2.37E-07 2.38E-07 4.46 4.17 4.01
DG- 20 1.54E-04 2.18E-05 2.20E-05 5
30 1.79E-05 2.46E-06 2.13E-06 5.32 5.37 5.75
40 3.79E-06 5.35E-07 5.18E-07 5.39 5.31 4.92
50 1.11E-06 1.61E-07 1.46E-07 5.50 5.39 5.69
DG- 10 9.72E-04 1.59E-04 2.00E-04 6
20 1.56E-05 2.13E-06 2.14E-06 5.96 6.22 6.55
30 1.14E-06 1.64E-07 1.91E-07 6.45 6.33 5.96
40 2.17E-07 2.97E-08 3.59E-08 5.77 5.93 5.82
DG- 5 2.24E-02 4.15E-03 3.11E-03 7
10 1.76E-04 2.75E-05 2.86E-05 6.99 7.24 6.76
20 1.67E-06 2.28E-07 2.26E-07 6.72 6.91 6.98
25 3.60E-07 4.96E-08 6.27E-08 6.86 6.84 5.74
DG- 5 5.50E-03 1.22E-03 1.46E-03 8
10 4.63E-05 6.26E-06 6.95E-06 6.89 7.61 7.71
15 1.62E-06 2.20E-07 2.29E-07 8.28 8.26 8.42
20 2.05E-07 2.80E-08 2.28E-08 7.18 7.17 8.01
DG- 4 9.11E-03 1.80E-03 3.44E-03 9
8 4.97E-05 7.51E-06 6.93E-06 7.52 7.90 8.96
10 7.50E-06 1.05E-06 1.18E-06 8.47 8.81 7.95
15 2.40E-07 3.34E-08 3.09E-08 8.49 8.51 8.98
DG- 4 3.95E-03 7.89E-04 1.42E-03 10
8 1.01E-05 1.44E-06 1.52E-06 8.61 9.09 9.87
10 1.44E-06 2.00E-07 2.27E-07 8.74 8.85 8.51
12 2.67E-07 3.70E-08 3.77E-08 9.26 9.25 9.85
Table 2: and errors and convergence rates for the 2D isentropic vortex problem for the ADER-DG- scheme supplemented with a posteriori ADER-WENO3 subcell limiter for to from top to bottom.

5.3 Shu-Osher oscillatory shock tube

This test shuosher2 is a 1D hydrodynamic shock tube. The downstream flow has a sinusoidal density fluctuation with a wave length of and an amplitude of . A Mach shock front is initially located at on domain . The left and the right states are given by , , and , and . The final time is set to . This problem involves small scales after the shock has interacted with the sine wave that can be captured either with a fine enough mesh or with high order accurate method. Here a very coarse mesh made of cells in -direction and in -direction is chosen. The numerical flux was the Osher-type flux presented in OsherUniversal .
The results of ADER-DG- with a posteriori WENO3 subcell limiter method are presented in Figure 5. A 1D view is first depicted, where each polynomial is represented by sample points per cell. A reference solution (straight line) obtained by a classical third order ADER-WENO scheme on a very fine mesh is also plotted. The quality of the result is excellent for a cell mesh. The bottom part of the figure displays the cell centered polynomials. The red cells have been limited, therefore updated with WENO3 on subcells whereas the blue cells have been updated with unlimited ADER-DG- on the main grid. It is important to note that the DG polynomials are almost continuous across all cell boundaries, except for the shock wave, where a real discontinuity is preserved. In other words the unavoidable limiting has not smeared the shock wave on two cells of the main grid, but it has been able to keep the subcell description of this discontinuity, thus leading to a shock wave spread over only one or two subcells. Note that the leftmost waves are actually contained within one single cell. These waves are steepening and will later become genuine shock waves222This is the same effect for Burgers equation when a smooth initial profile evolves into a shock wave.. As expected the a posteriori detection procedure has flagged these waves for subsequent subcell correction, hence their red color.

Figure 5: Shu-Osher problem at with ADER-DG- with a posteriori ADER-WENO3 subcell limiter on a mesh — Top: 1D cut through the numerical solution (symbols) vs reference solution (ultra fine ADER-WENO solution in straight line). Any polynomial is represented by sample points per cell — Bottom: troubled cells, which have been updated with ADER-WENO3 on the subgrid, are shown in red, blue cells have used the unlimited ADER-DG- scheme on the main grid.

5.4 Double Mach reflection problem

Next we have run the 2D double Mach reflection problem of a strong shock that was proposed in woodwardcol84 . This test problem involves a Mach shock in a perfect gas with which hits a ramp with the -axis. Using Rankine-Hugoniot conditions we can deduce the initial conditions in front of and after the shock wave

(33)

where is the coordinate in a rotated coordinate system. Reflecting wall boundary conditions are prescribed on the bottom and the exact solution of an isolated moving oblique shock wave with shock Mach number is imposed on the upper boundary. A Rusanov (local Lax-Friedrichs) flux has been used for this test problem. Inflow and outflow boundary conditions are prescribed on the left side and the right side, respectively.
The computational domain is given by and the main grid is built using a characteristic length of , leading to computational cells. We solve this problem with four schemes: ADER-DG- with , , and , all supplemented with the a posteriori ADER-WENO3 subcell limiter. The results are depicted using isolines ranging from to for the density variable at , see Fig 6. We first present the density contour lines along with the cells colored by the indicator function : red cells have been recomputed using the subcell ADER-WENO3 scheme (), whereas blue cells () have kept the original unlimited ADER-DG- method. The shock waves are well detected for all DG schemes. Moreover, the interaction zone where a lot of small scale vortex structures develop is almost free of limiting for ADER-DG- and ADER-DG- and more troubled cells are detected by our a posteriori method for a low order ADER-DG- scheme, which is quite remarkable. One would have intuitively expected that higher order schemes tend to create more oscillations and thus need more limiting. It seems, however, that the higher order DG scheme has better subcell resolution capabilities and therefore can treat these vortex structures as smooth subscale features, without the need of a limiter. Also more erroneously detected troubled cells are observed for the ADER-DG- scheme in the upper part of the figure, behind the main shock wave. Only minor details seem to differ in the indicator function between ADER-DG- and ADER-DG-. Nonetheless, the numerical structures represented by the density isolines seem to be richer with the ADER-DG- scheme.

Overall, also for this very difficult problem which exhibits at the same time strong shock waves and smooth flow features, our detection procedure seems to identify discontinuities associated with shock waves properly and it ignores complex but smooth vortex-type flow structures, as expected.

Figure 6: Double Mach reflection problem at with ADER-DG-, ADER-DG- and ADER-DG- (from top to bottom) supplemented with a posteriori ADER-WENO3 subcell limiter — isolines ranging from to are displayed for the density variable (black lines) along with the troubled cells (red) and the unlimited cells (blue) at this time level. The characteristic mesh spacing is only in all cases.

In order to confirm this observation we propose in Fig 7 two zooms on the interaction zone for ADER-DG- schemes with , , , and . The main grid is also plotted in light gray color to compare the numerical thickness of a shocks wave and the size of the vortex-type structures with respect to the main grid size. From this figure we confirm that the DG scheme behaves better with increasing , the degree of the DG polynomial. ADER-DG- dissipates less and, as such, is capturing more structures than the other schemes. Also the smearing of the shocks is less pronounced with ADER-DG- than with ADER-DG- on the same grid. All schemes capture the shock waves properly without spurious oscillations, but the numerical dissipation of the lower order schemes destroys the small scale vortex structures that are obtained on this rather coarse grid only by the highest order DG schemes. The limiter does not seem to destroy this ability. In other words, the subcell resolution property of the DG scheme seems to be very well preserved.

Figure 7: Double Mach reflection problem at with ADER-DG-, , and (from top left to bottom right) supplemented with a posteriori ADER-WENO3 subcell limiter — Zooms on the interaction zone. isolines ranging from to are displayed for the density variable (black lines). The mesh is underlined in gray color.

From the results of this test case we can conclude that the a posteriori subcell limiter provides a valid detection strategy of troubled cells and that the ADER-WENO3 scheme used on the subgrid is able to preserve the overall accuracy of the high order DG scheme on the main grid.

5.5 Forward facing step

In this section we simulate the so called forward facing step (FFS) problem, also proposed by Woodward and Colella in woodwardcol84 . It consists of a Mach 3 wind tunnel with a step. The initial condition consists in a uniform gas with density , pressure , velocity components and . The computational domain is given by and reflective wall boundary conditions are applied on the upper and lower boundary of the domain, whereas inflow and outflow boundary conditions are applied at the entrance and the exit. The solution of this problem involves shock waves interacting with the boundaries and it is run up to a final time of . The mesh is uniform, using a characteristic length of , leading to about cells on the main grid. For this test we run the ADER-DG- scheme with a Rusanov flux, supplemented with the a posteriori ADER-WENO3 subcell limiter (SCL). In Fig. 8 we display the density variable at the final time. The top panel represents an extruded density (the azimuth and the color represent the numerical density). From this view we can clearly see several smooth zones separated by the shock waves. Furthermore, we can distinguish the Kelvin-Helmholtz instabilities developing along the top shear wave. These instabilities are maintained when they move across the shock waves. Moreover, the unsteady vortices create small amplitude acoustic waves, which are expected from such unsteady flows. Yet, they are often smeared by more dissipative numerical schemes. On the same figure we propose in the middle panel the classical black and white isolines in contrast with the previous figure. It is clear that the perturbed isolines emanating from the vortices are due to the acoustic waves, previously seen on the extruded 3D view. Finally the bottom panel shows the limited cells updated with the subcell WENO3 scheme (in red) and the unlimited ADER-DG- cells (in blue). As expected, the great majority of the computational domain is computed with the unlimited ADER-DG- scheme, apart from the shock waves, where the subcell limiter mechanism seems to act properly.

Figure 8: Forward Facing Step problem at with ADER-DG- supplemented with a posteriori ADER-WENO3 subcell limiter — Top-panel: extruded 3D density (color and azimuth) — Middle-panel: density represented as isolines — Bottom-panel: troubled cells (red) updated using the subcell ADER-WENO3 method and unlimited ADER-DG- cells (blue).

5.6 2D Riemann problem

In this section we consider a set of two-dimensional Riemann problems which has been proposed and extensively studied in schulzrinne ; kurganovtadmor and which has also been recently used to construct genuinely multidimensional Riemann solvers of the HLL type, see balsarahlle2d ; balsarahllc2d ; BalsaraMultiDRS . The computational domain is and the initial conditions are given by

(34)

The initial conditions and the final simulation time, , for the five configurations presented in this article are listed in Table 3. For more information about the other configurations the reader is referred to schulzrinne ; kurganovtadmor . For the simulation we have used a main grid of elements with a DG- scheme supplemented with our a posteriori WENO3 subcell limiter.

#
RP1 0.5323 1.206 0.0 0.3 1.5 0.0 0.0 1.5 0.25
0.138 1.206 1.206 0.029 0.5323 0.0 1.206 0.3
RP2 0.5065 0.8939 0.0 0.35 1.1 0.0 0.0 1.1 0.25
1.1 0.8939 0.8939 1.1 0.5065 0.0 0.8939 0.35
RP3 2.0 0.75 0.5 1.0 1.0 0.75 -0.5 1.0 0.30
1.0 -0.75 0.5 1.0 3.0 -0.75 -0.5 1.0
RP4 1.0 -0.6259 0.1 1.0 0.5197 0.1 0.1 0.4 0.25
0.8 0.1 0.1 1.0 1.0 0.1 -0.6259 1.0
RP5 1.0 0.7276 0.0 1.0 0.5313 0.0 0.0 0.4 0.25
0.8 0.0 0.0 1.0 1.0 0.0 0.7276 1.0
Table 3: Initial conditions for the 2D Riemann problems numbered from to . These further correspond to Configurations 3, 4, 6, 8 and 12 in kurganovtadmor

In Figs. 9-10 the numerical results are presented for the first three and for the last two Riemann problems, respectively. In the left panels we show the distribution of the density at the final time, with equidistant isolines between the minimum and the maximum value. The number of isolines is for RP1 and RP2, for RP3 and RP4, and for RP5. In the right panels we show, as usual, the corresponding mesh and the cells, colored in red, which have been updated with the subcell WENO3 scheme, while the unlimited cells are marked in blue. It seems that the limiter is active only along strong discontinuities. Apart from these waves (and some parasitical cells along with the start-up error for RP3) the limiter is inactive in smooth regions, leading to an optimal precision given by the unlimited ADER-DG- scheme.
The computational results, and in particular the generation of the main structures in all these two dimensional Riemann problems, are in good agreement with the literature kurganovtadmor . Nonetheless, in these simulations the numerical dissipation is drastically reduced and, as such, the numerical solution shows much more small-scale features than usually reported, which are for example due to the Kelvin-Helmholtz instability along shear waves (see RP3 and RP4 in Figs. 9-10 respectively).

To have a reliable comparison, we have repeated the test case RP3 using a sixth order finite volume ADER-WENO scheme with space-time adaptive mesh refinement (AMR), using the method described in AMR3DCL . The results of this comparison are illustrated in Fig. 11: in the left panel we have simply reported again the solution of Fig. 9, while in the right panel we can see the solution with the sixth order ADER-WENO scheme using AMR. In the AMR simulations, the level 0 grid is composed of elements which has then been adaptively refined using two levels of refinement and a refinement factor of 5 (right panel), which leads to an equivalent resolution of elements on the finest level. The high order finite volume simulation with an AMR code clearly confirms the onset of a Kelvin-Helmholtz instability, which means that the results obtained with the high order ADER-DG scheme with subcell limiter are reliable, since this kind of physical instability actually should be observed in numerical simulations for sufficiently fine meshes or sufficiently high order accurate schemes.

For these classical test cases it seems that the coupling of a posteriori subcell limiter and high order DG schemes is a valid tool to capture the discontinuous waves without spurious oscillations, and also the smooth part of the flow without excessive numerical dissipation.

Figure 9: 2D Riemann problems simulated with an ADER-DG- method supplemented with a posteriori ADER-WENO3 subcell limiter — Left panels: isolines of the density variable — Right panels: limited cells (red) and unlimited cells (blue). From top to bottom RP1, RP2 and RP3.
Figure 10: Same caption as Fig.9, but for RP4 and RP5 (from top to bottom).
Figure 11: Comparison of RP3 solved with two different numerical methods. Left panel: same as Fig.9, obtained with the ADER-DG- method with subcell limiter on a coarse main grid. Right panel: solution obtained with a sixth order ADER-WENO6 finite volume scheme with AMR (equivalent resolution: ), showing the Kelvin-Helmholtz instability on the shear waves.

5.7 Shock-vortex interaction

As a final test in two space dimensions, we have considered the interaction of a vortex with a steady shock wave. Originally proposed by donat , this test is a true benchmark for a high order numerical scheme, as it involves a complex flow pattern with both smooth features and discontinuous waves. The initial conditions, defined over the computational domain , are given by a stationary normal shock wave placed at and by a vortex, which is placed at . The shock Mach number is denoted by and inside the vortex we have the following distribution of the angular velocity: