1. Introduction
Highorder numerical methods have been successfully applied to a wide range of applications [45, 47, 50]
. These methods promise higher accuracy per degree of freedom when compared with traditional loworder methods, and have the potential to achieve high efficiency on modern computing architectures
[18, 3, 9]. For example, in the field of computational fluid dynamics, such methods have seen particular success when applied to underresolved turbulent flows, such as in the context of implicit large eddy simulation (ILES) [34, 38, 4]. However, a critical issue that must be addressed when applying highorder methods to convectiondominated problems is their robustness, especially in the context of nonlinear problems with discontinuous features such as shock waves [35, 52, 37].In particular, discontinuous Galerkin (DG) methods have seen considerable success when applied to convection dominated problems [7]. These methods possess many desirable properties, such as arbitrary formal order of accuracy, and suitability for use with unstructured meshes. The robustness of DG methods is the subject of a large body of research [22, 34]. For scalar conservation laws and symmetric systems, the DG method satisfies a cell entropy inequality [19]. However, for general hyperbolic systems, such as the Euler equations, techniques such as flux differencing are required to ensure entropy stability [8, 6, 5, 39]. Furthermore, the use of high degree polynomials can introduce oscillations, and therefore limiters or artificial viscosity techniques are often used for bounds preservation, monotonicity, and shock capturing [43, 24, 25, 42].
An alternative to the above stabilization and limiting strategies is an approach developed by Guermond, Popov, and colleagues, based on invariant domain preserving (IDP) discretizations and convex limiting [15, 13, 16, 27]. A desirable property for numerical discretizations of hyperbolic conservation laws is invariant domain preservation: if the exact solution to the conservation law lies in a convex invariant set, then the numerical solution should as well [15]. This is a generalization of the concept of a discrete maximum principle, and will ensure that the the discretization is bounds preserving, positivity preserving, and nonoscillatory. Suitable loworder invariant domain preserving (IDP) discretizations have been paired with highorder discretizations using convex limiting or algebraic flux correction strategies to obtain secondorder accurate methods that preserve specific invariant domains [13, 16, 27].
In this work, we develop highorder discontinuous Galerkin methods that satisfy invariant domain preserving properties using a convex limiting strategy. The limiting strategy makes use of a novel sparse loworder IDP method whose stencil does not grow with the polynomial degree of the corresponding highorder method. Crucially, the accuracy of the loworder method does not degrade as the polynomial degree of the highorder method is increased, as is observed to occur with more naive graph visocity approaches. Related strategies for sparsifying the convective operator for Bernstein basis finite element methods were previously developed by Kuzmin and colleagues [28]. The fluxcorrected method is obtained by performing an efficient, dimensionbydimension subcell flux correction procedure, blending the loworder IDP method with the highorder target DG method. The resulting method is conservative, and can satisfy any number of constraints on quasiconcave functionals specified by the user (cf. [13]). Since the accuracy of the loworder method does not decrease with the polynomial degree of the target method, we observe more accurate results using higherorder methods with a fixed number of degrees of freedom, even on problems with discontinuous solutions. This method can also be combined with a subcell resolution smoothness indicator to alleviate peak clipping effects near smooth extrema.
The structure of this paper is as follows. In Section 2, we formulate the highorder DG discretization, and state some key properties. In Section 3, we introduce a new loworder sparsified discretization that can be rendered invariant domain preserving using a graph viscosity approach. We develop a subcell flux correction strategy for blending the highorder (target) method and the loworder IDP method in Section 4. As specific examples, applications to the linear advection equation with variable velocity field and the Euler equations of gas dynamics are discussed. A number of numerical test cases demonstrating the effectiveness of the method on both scalar equations and systems of hyperbolic conservation laws are presented in Section 5. Finally, we end with some concluding remarks in Section 6.
2. Governing equations and discretization
Consider a system of hyperbolic conservation laws,
(1) 
with solution , . The flux function is given by . The spatial dimension is denoted , and the number of solution components is . The initial conditions are given by . Closely associated with the problem (1) is the following one dimensional Riemann problem, which will be important both for the definition of invariant domain preservation, and for the formulation of the discontinuous Galerkin discretization. Let be a pair of admissible states, and let
by any unit vector. We assume that the Riemann problem
(2) 
has a unique selfsimilar entropy solution . We denote by the maximum wave speed for (2), for which we have if , and if .
Definition 1.
For example, the maximum principle implies that any interval is an invariant set for scalar conservation laws. For the Euler equations, the set of states with positive density, positive internal energy, and satisfying a minimum principle on specific entropy is a convex invariant set. Additional examples of invariant sets for systems of conversation laws are given in [15, 16]. It will be desirable to construct discretizations of (1) that are invariant domain preserving (IDP), meaning that if the approximate solution lies in a convex invariant set at some time , then the solution will remain in for all time .
The strategy we present here for developing IDP discretizations for (1) is as follows. We first formulate a highorder discontinuous Galerkin (DG) discretization that will serve as a target scheme. This discretization will in general not be invariant domain preserving. We will then modify this highorder discretization to generate a robust loworder discretization. These modifications take the form of first sparsifying the method to reduce the size of the stencil, and then adding a graph viscosity term (cf. [15]), which guarantees that the resulting discretization is invariant domain preserving. Finally, a subcell flux corrected transport (FCT) technique is used to blend the loworder IDP method and the highorder target method in such a way that specified convex invariant sets are preserved.
2.1. DG formulation
We begin by defining the highorder DG discretization for equation (1). The spatial domain
is discretized with a mesh of tensorproduct elements denoted
. Each element is the image of the reference element (the unit cube in dimensions) under a transformation mapping . To define the standard discontinuous Galerkin finite element space , first consider the space defined on the reference element consisting of all multivariate polynomials of degree at most in each variable. On a given element , we define the space to be spanned by functions , where is the element mapping, for all . Then, the space is defined as(3) 
Note that no continuity is enforced between adjacent elements. To represent approximate solutions to (1), we also consider the vector version of this space .
We use a nodal GaussLobatto basis for the space . Let denote the GaussLobatto points in the interval , and let
denote the Lagrange interpolating polynomial satisfying
, where is the Kronecker delta. These functions form a basis for the space in one dimension. The basis for is formed by taking the tensor product of the onedimensional basis. To be precise, we define a function , where denotes the multiindex , and . This basis can also be seen to be the nodal interpolation basis corresponding to the Cartesian product of the onedimensional GaussLobatto nodes.We approximate the solution to (1) by , multiply the equation by a test function , and integrate over the domain , integrating the flux term by parts over each element . Because the space is discontinuous, the fluxes are not welldefined on element interfaces. Consider two neighboring elements, and . Let denote the trace of from within , and similarly for . We therefore introduce a singlevalued numerical flux function , obtaining the weak formulation
(WF) 
Integrating the second term on the lefthand side once more by parts, elementbyelement one obtains what is know as the strong formulation,
(SF) 
Note that at the continuous level, the formulations (WF) and (SF) are equivalent. However, after discretization, they may differ because of inexact integration.
For the purposes of discretization, it is convenient to transform the integrals in both (WF) and (SF) to integrals over the reference element . This is done using a standard transformation of the governing equation (1) from each element to the reference element . Consider a given element with transformating mapping . Let denote the Jacobian matrix of the mapping . The inverse of the Jacobian is used to define the contravariant fluxes
(4) 
Then, on the reference element, the solution evolves according to the transformed conservation law
(5) 
where .
2.1.1. Numerical flux functions
An important aspect of DG methods is the choice of numerical flux function . The numerical flux functions are typically chosen to be either exact or approximate Riemann solvers for the onedimensional Riemann problem in the normal direction at element interfaces. In this work, we will make use of the simple local LaxFriedrichs numerical flux function. The reason for this choice is that the LaxFriedrichs flux is compatible with the graph viscosity used to ensure that the low order discretization is invariant domain preserving, as will be discussed in greater detail in Section 3. The LaxFriedrichs flux is defined by
(6) 
where is an upper bound for the maximum wave speed of the Riemann problem
(7) 
2.2. Collocated and onedimensional operators
The discontinuous Galerkin spectral element method (DGSEM) is distinguished from other DG methods by a specific choice of quadrature rule and basis. We proceed by choosing a nodal basis for the space and then approximating the integrals in (WF) and (SF) with collocated quadrature rules. Typically, GaussLegendre or GaussLobatto nodes are chosen for the basis functions. In this work, we will solely make use of GaussLobatto nodes and quadrature.
Due to the use of tensorproduct basis and quadrature, the DG operators possess a Kroneckerproduct structure [40]. Since the nodal points and quadrature points are collocated, there is no need for an interpolation operator. The onedimensional mass matrix on the reference element, denoted , is given by a diagonal matrix with quadrature weights on the diagonal. The weighted onedimensional differentiation matrix is obtained by evaluating the derivatives of the basis functions at the nodal points and multiplying by quadrature weights,
(8) 
We will often make use of the following two simple properties of the differentiation matrix:
Proposition 1.
The weighted differentiation matrix satisfies the following two useful properties, which we will make use of extensively in this work:
(P1)  
(P2) 
On the reference element, the local mass and differentiation operators can be obtained through Kronecker products. For instance, for we have
(9) 
We can also define the left endpoint evaluation matrix , which is zero except for the first entry of the diagonal, which takes value one, and likewise the right endpoint evaluation matrix , which is zero except for the last entry of the diagonal. Given these definitions, integrals over the boundary of the reference element can be computed using Kronecker products,
(10)  
2.3. Metric terms and transformed operators
Having defined the mass, differentiation, and boundary operators on the reference element as above, we now wish to transform these operators to act in physical coordinates. Consider a fixed element with transformation mapping and Jacobian . Let denote the values of the Jacobian matrix of evaluated at the GaussLobatto nodal points. These terms are evaluated using the freestreampreserving procedure descried in [23]. Let denote the diagonal matrix whose entries are given by . Similarly, let denote the ( block) diagonal matrix whose entries are given by . Notice that the elemental mass matrix corresponding to the element is given by .
Let denote the vector of coefficients (i.e. nodal values) of on and let denote the vector of values of evaluated at the nodal points (i.e. ). For each face , let denote the values of evaluated at the dimensional GaussLobatto nodes on face (where the trace of is taken from within ), dotted with the scaled normal vector facing outwards from . Likewise denotes the nodal values of .
2.4. Conservation and constant preservation
The governing equation (1) satisfies the following two simple properties:

(Conservation). Assuming periodic or compactly supported boundary conditions, .

(Constant preservation). If is spatially constant, then . In particular, if depends only on (and not, for example, on the spatial variable ), then if being spatially constant implies that .
We would like the discretization to satisfy the analogous properties at the discrete level.
2.4.1. Conservation
First, we consider conservation. The analogous statement at the discrete level is that
(16) 
where is a vector of all ones. First, note that since the numerical flux function is singlevalued, we have
(17) 
The remaining terms in both the strong form and weak form are local to each element and do not involve contributions from face neighbors. In the weak form (14), we see that by property (P1). We can then conclude that
(18) 
proving the conservation property for the discretized weak form.
For the discretized strong form, we make use of property (P2). For any
(19) 
and we therefore obtain conservation for the strong form as well:
(20) 
We summarize the above arguments in the following proposition:
2.4.2. Constant preservation
Now we turn to the property of constant preservation. We supposed that the flux is spatially constant, and therefore is also everywhere equal to the same constant. As a result, this property is easily proven for the strong form of the discretization: , and so the third term on the lefthand side of (15) is zero. Therefore, the discretization will preserve constants if the following identity holds:
(22) 
This socalled metric identity is discussed in detail in [23], and can be enforced through proper evaluation of the entries of . In this case, we have the following result:
Proposition 3.
Let satisfy the discretized strong form (15). Furthermore, suppose that is spatially constant. Then, .
3. Construction of the IDP loworder method
We now modify the discretization (15) with the goal of obtaining a method which is invariant domain preserving (IDP). In Guermond and Popov [15] this is done by adding a graph viscosity term that is based on the guaranteed maximum speed (GMS) of the hyperbolic system. One potential issue with this approach is that the amount of graph viscosity added to the discretization increases as the size of the discrete stencil increases. As a result, applying this technique to highorder methods with large stencils results in very dissipative methods and typically poorquality results, as observed in [33], and further illustrated in Section 3.4. In order to address this issue, we are interested in creating an IDP discretization that is compatible in a certain sense with (15), yet based on a more compact stencil. Because the graph viscosity term is itself firstorder accurate, the underlying sparse discretization is not required to be highorder accurate.
We make the following simple modification to the DGSEM method described above. Notice that the wide stencil of the highorder method is a consequence of the fact that the onedimensional differentiation matrix is dense. We therefore replace with a sparser version that is firstorder accurate. is obtained by integrating the derivatives of piecewise linear basis functions on the mesh defined by the GaussLobatto points in the interval . is therefore given by
(23) 
Then, for each , we construct operators by replacing with in the definition (9). It is easy to see from (23) that also satisfies the row and columnsum properties (P1) and (P2). The modified operators then replace the standard DGSEM derivative operators in the formulation (14), in order to obtain the following modified discretization:
(24) 
3.1. Conservative correction
For our purposes, it is important to ensure that the modified formulation (24) satisfies the conservation and constant preservation properties described in section 2.4. Because (14) is based on the weak formulation, the conservation property follows immediately from the zerosum property (P1) (and the fact that the numerical fluxes are unique). However, the constant preservation property will not hold in general. The reason for this is that the modified differentiation matrix cannot differentiate exactly the metric terms in , and therefore the modified discrete metric identities no longer hold.
To remedy this issue, we make use of modified metric terms which are perturbations of the highorder metric terms , but are designed such that the firstorder discretization satisfies the metric identities. Since the modified loworder method is itself only firstorder accurate, using an perturbation of the metric terms is acceptable.
The modified metric terms are constructed as follows. We set , where are block diagonal correction matrices. For each element, we would like to enforce the identity
(25) 
where denotes the th component of the outwardfacing normal at the GaussLobatto points of face . This is equivalent to the following underdetermined system of equations:
(26) 
Solving for the correction matrices , we obtain the system
(27) 
Note that the differentiation matrix on the lefthand side has a nullspace consisting of all constant functions. In order for a solution to exist, the righthand side must be orthogonal to this nullspace. In other words, the entries of the righthand side vector must sum to zero. This can be seen to be true for the second term on the righthand side by a simple application of the zerosum property (P1). In order to show this property for the first term on the righthand side, we make use of the summationbyparts property (19):
(28) 
Then, the metric identity (22) implies that this term is equal to zero. Therefore a solution to (27) exists. On each element, the minimum norm solution to (27
) is computed for the perturbed metric terms. This can be performed using e.g. the singular value decomposition, which is performed once as an offline precomputation on the reference element, and then simply reused for each element
. The righthand side (as well as the metric terms themselves) scale as , and therefore the entrywise error satisfies .Using the above procedure to construct the modified metric terms, we obtain the loworder discretization
(29) 
As a consequence of the choice of metric terms, this discretization satisfies the following properties.
Proposition 4.
Let satisfy (29). Then, the following two properties hold:

Conservation:

Constant preserving: if is spatially constant then .
Next, the discretization (24) is modified to render the resulting method invariant domain preserving (IDP). Because of the local nature of the DG method, and because of the choice of LaxFriedrichs numerical flux (cf. Section 2.1.1), it is possible to perform this modification in an entirely local fashion. First, we rewrite (24) in a form that will be more convenient for our purposes. Note that by definition of the LaxFriedrichs flux (equation (6)), the term of the form can be written as
(30) 
where is the number of faces on which the th node lies (i.e. depending on if is an interior node or if it lies on a face, edge, or corner of the element), and are weighted normal vectors, and is the index set consisting of all nodes that are face neighbors of . The vector is equal to the normal vector evaluated at node , pointing outwards from the element to which node belongs, weighted by the element of surface area and face quadrature weight. If the node does not lie on an element face, then and , and in this case we leave undefined. Inserting this expression into (24) and writing the volume terms in terms of coefficients , we obtain the following equivalent formulation:
(31) 
where is the th diagonal entry of the mass matrix, and is the set of all indices in the stencil of within the same element. For any , we define the symmetric graph viscosity coefficients by
(32) 
where . For convenience of notation, we will use the convention that . Having defined the quantities , the discretization (31) is modified with the addition of a viscosity term
(33) 
The above expression is simplified by combining the volume and boundary terms. Let , and let , where denotes the corresponding coefficient of from the boundary terms. Here we use the convention that if and likewise if . At this point, notice that the LaxFriedrichs term may be combined with the graph viscosity by appropriately defining the viscosity coefficients. Let for be given by
(34) 
Note that by this definition is equivalent to replacing with in (32). Using these definitions, (33) simplifies to
(35) 
Proposition 5.
The discretization given by (35) satisfies the following properties.

Since the loworder discretization is constant preserving,

By symmetry of the coefficients , the graph viscosity contributions sum to zero:
3.2. Invariant domain preservation
We now set out to prove that the discretization defined by (35) is invariant domain preserving (IDP). We will make use of a strong stability preserving (SSP) RungeKutta method for the temporal discretization [11]. Such methods can be written as a convex combination of forward Euler steps, and so it suffices to prove the IDP property for a forward Euler update. We therefore consider the forward Euler discretization of (35), given by
(36) 
As in [15], we use the property to rewrite (35) in the form
(37) 
At this point we can introduce the socalled “bar states” (also referred to as “intermediate limiting states”, cf. [15]), which are defined by
(38) 
Proposition 6.
Suppose is a convex invariant set of (1) such that Then, the bar states belong to .
Proof.
See [16]. ∎
Having defined the bar states, we rewrite (37) as a convex combination
(39) 
As an immediate consequence of this rewriting is the following proposition.
Proposition 7.
Suppose the following CFL condition holds:
(40) 
where . Let denote a convex invariant set of (1), such that for all . Then, for all .
3.3. Application: linear advection
Consider the linear, scalar advection equation
(41) 
where is a prescribed velocity field. We assume that the velocity field is divergence free, i.e. . We are interested in ensuring that the loworder method is bounds preserving, i.e. if for all , then for all .
If is spatially constant, then the method described above applies immediately to this case, and the loworder method defined by (36) is bounds perserving. However, if is spatially variable, then some modifications to the above method are required. The reason for this is that in this case, Proposition 6 may fail to hold. For instance, if is spatially constant (i.e. for all , for some fixed ), then we would expect to be equal to the same constant. However, since is spatially varying, we have, in general, and therefore .
To avoid this issue, we slightly modify the formulation, using an approach similar to that developed by Kuzmin in [27]. We define the modified bar states for the advection equation by
(42) 
where , and denotes the coordinates of the th node of the mesh. The modified graph viscosity coefficients for the advection equation are given by
(43) 
It is clear from this definition that if then . As a consequence, defining the modified update for as the convex combination of bar states given by (39), we see that , as desired.
Remark.
The property was essential to writing the loworder update in terms of the bar states. In the context of linear advection, the analogous property is . This is the discrete equivalent to the divergence free constraint, . In order to ensure that this property holds, the conservative correction procedure described in Section 3.1 must be modified to take into account the velocity field .
3.4. Comparison with unsparsified method
As mentioned in Section 3, the motivation for introducing the sparsified derivative operators is that the addition of the graph viscosity term to the “unsparsified” operators can result in overly dissipative results when the stencil size is increased.
To illustrate this point, we consider the advection of a sine wave in one dimension. The domain is taken to be , the initial condition is , and periodic boundary conditions are enforced. We integrate the equation until a final time of , at which point the solution is identical to the initial condition. We compare the sparsified loworder IDP method to the graph viscosity method applied to the unsparisified DGSEM operator. We fix the total number of degrees of freedom to be 128, and consider polynomial degrees 3, 7, 15, and 31 (corresponding to mesh sizes of 32, 16, 8, and 4 elements, respectively). The results are shown in Figure 1. It is immediately clear that increasing the polynomial degree causes a large degradation in quality of the unsparsified methods. This is because the size of the stencil grows with the polynomial degree, causing each degree of freedom to be coupled to other degrees of freedom. Each such connection results in an additional graph viscosity contribution, rendering the method overly diffusive. On the other hand, the sparsified operator has a stencil size of , and as a consequence, we do not observe a degradation of the results with increased .
4. Flux correction and limiting strategies
We now combine the loworder IDP discretization (35) with the highorder discretization (15), to obtain a boundpreserving scheme. To this end, we consider a forward Euler time discretization of both the loworder and highorder approximations (as mentioned above, the extension to highorder time integration is straightforward using SSP methods). The provisional highorder update for element is given by
(44) 
and the loworder IDP update is given by
(45) 
We introduce the shorthand notation
(46) 
Note that the provisional highorder solution may not satisfy convex invariants, such as maximum principles, positivity, etc. However, the loworder solution is guaranteed to be IDP as long as the CFL condition (40) is satisfied. We then define the fluxcorrected update by
(47) 
where is a limiting factor that is yet to be determined.
To ensure conservation of the flux corrected solution, it is important that the correction terms do not impair the following zerosum property:
Proposition 8.
Let denote the set of indices associated with the element . Then,
(48) 
Proof.
In fact, the above proposition can be extended to consider “lines” of degrees of freedom in a dimensionbydimension fashion. Recall that each element contains nodes, which we can index as . Fix a dimension . Then, let denote the set of all indices of nodes in the element , with local index , where the th index of has been replaced by with . Furthermore, the residuals within a given element are decomposed into contributions corresponding to each coordinate dimension:
(49) 
Proposition 9.
Consider a “line” of nodes within an element (for any and multiindex as described above). Then
(50) 
Proof.
This property follows from the Kronecker product definition of the operators and , e.g. as written in (9). ∎
4.1. Linear constraints for scalar problems
Suppose the governing equation is a scalar conservation law. In this case, the equation satisfies a local maximum principle, and any interval is a convex invariant set. Since the loworder method is IDP, we have that if for all , then . We wish to enforce a similar property for the highorder limited quantity . Specifically, for each , we choose bounds and and enforce . Since this constraint is linear (as opposed to nonlinear constraints such as entropy inequalities), the limiting procedure is relatively simple.
4.1.1. Bounds
In order to choose the limiting factors , we must first choose bounds and for all . We can consider several choices of bounds:

,

,

.
We have shown earlier in this document that and preserve convex invariant sets of the governing equation. Therefore, these quantities satisfy the discrete maximum principle. Note that is a convex combination of , and hence represents a more restrictive bound. To define and we may take the minimum and maximum of some combination of these quantities. For the remainder of this section, we choose the bounds naturally satisfied by the loworder discretization:
(51) 
4.1.2. Elementwise limiting for linear constraints
A Zalesaktype limiter can be used to compute a provisional limiting factor for each , cf. [33, 53]. Then, can be defined as the minimum of all over the element :
(52) 
Given this definition, is constant on each element, and therefore Proposition 8 implies that the resulting FCT method is conservative. However, choosing to be constant over an entire element may be overly pessimistic if the polynomial degree is high, and therefore we consider also subcell limiting.
4.1.3. Subcell limiting for linear constraints
In order to improve the resolution of the fluxcorrected solution, we consider subcell limiting, where the correction factors are allowed to vary within each element. The idea of the subcell limiting is closely related to that presented in [33, Section 4.5]. In that work, a onedimensional subcell limiting algorithm was described, and the extension to multiple dimensions was proposed as a minimization problem. Instead of solving a minimization problem, Proposition 9 allows for the use of the onedimensional subcell limiting procedure along lines of nodes, in a dimensionbydimension fashion. For simplicity of notation, let denote the antidiffusive flux at the th node in the th dimension: .
We consider a decomposition of a given element into subcells, as illustrated in Figure 2. Note that the subcell boundaries are placed in between nodes, such that every node can be considered as a “subcellcentered value”. We consider the set of all “subcell faces,” which are the set of all interior subcell faces (the thin blue lines in the figure). For each subcell face, we assign an antidiffusive flux that is obtained by summing the nodal antidiffusive fluxes lying on one side of the face.
To introduce notation, we refer to the subcell faces by a pair , where . The index indicates that the face is normal to the th unit vector. For a given face (indicated by a thick red line in Figure 2), let denote set of indices of nodes lying on one side of the face (indicated by red nodes in Figure 2). We then define the subcell face flux associated with the face by summing the directional nodal fluxes over the set :
(53) 
Fix a node and direction , and let and denote the subcell faces adjacent to in the th direction. Note that , and therefore, by Proposition 9, we can write the nodal antidiffusive flux as the difference of adjacent subcell fluxes:
(54) 
For each node , we introduce a nodal provisional limiting coefficient , obtained by limiting the sums of positive and negative parts of the adjacent subcell residuals,
Comments
There are no comments yet.