. These methods promise higher accuracy per degree of freedom when compared with traditional low-order 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 under-resolved 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 high-order methods to convection-dominated 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 . 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 . 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 . 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 non-oscillatory. Suitable low-order invariant domain preserving (IDP) discretizations have been paired with high-order discretizations using convex limiting or algebraic flux correction strategies to obtain second-order accurate methods that preserve specific invariant domains [13, 16, 27].
In this work, we develop high-order discontinuous Galerkin methods that satisfy invariant domain preserving properties using a convex limiting strategy. The limiting strategy makes use of a novel sparse low-order IDP method whose stencil does not grow with the polynomial degree of the corresponding high-order method. Crucially, the accuracy of the low-order method does not degrade as the polynomial degree of the high-order 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 . The flux-corrected method is obtained by performing an efficient, dimension-by-dimension subcell flux correction procedure, blending the low-order IDP method with the high-order target DG method. The resulting method is conservative, and can satisfy any number of constraints on quasiconcave functionals specified by the user (cf. ). Since the accuracy of the low-order method does not decrease with the polynomial degree of the target method, we observe more accurate results using higher-order 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 high-order DG discretization, and state some key properties. In Section 3, we introduce a new low-order sparsified discretization that can be rendered invariant domain preserving using a graph viscosity approach. We develop a subcell flux correction strategy for blending the high-order (target) method and the low-order 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,
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
has a unique self-similar entropy solution . We denote by the maximum wave speed for (2), for which we have if , and if .
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 high-order 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 high-order discretization to generate a robust low-order 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. ), which guarantees that the resulting discretization is invariant domain preserving. Finally, a subcell flux corrected transport (FCT) technique is used to blend the low-order IDP method and the high-order target method in such a way that specified convex invariant sets are preserved.
2.1. DG formulation
We begin by defining the high-order DG discretization for equation (1). The spatial domain
is discretized with a mesh of tensor-product 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
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 Gauss-Lobatto basis for the space . Let denote the Gauss-Lobatto 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 one-dimensional basis. To be precise, we define a function , where denotes the multi-index , and . This basis can also be seen to be the nodal interpolation basis corresponding to the Cartesian product of the one-dimensional Gauss-Lobatto 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 well-defined on element interfaces. Consider two neighboring elements, and . Let denote the trace of from within , and similarly for . We therefore introduce a single-valued numerical flux function , obtaining the weak formulation
Integrating the second term on the left-hand side once more by parts, element-by-element one obtains what is know as the strong formulation,
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
Then, on the reference element, the solution evolves according to the transformed conservation law
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 one-dimensional Riemann problem in the normal direction at element interfaces. In this work, we will make use of the simple local Lax-Friedrichs numerical flux function. The reason for this choice is that the Lax-Friedrichs 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 Lax-Friedrichs flux is defined by
where is an upper bound for the maximum wave speed of the Riemann problem
2.2. Collocated and one-dimensional operators
The discontinuous Galerkin spectral element method (DG-SEM) 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, Gauss-Legendre or Gauss-Lobatto nodes are chosen for the basis functions. In this work, we will solely make use of Gauss-Lobatto nodes and quadrature.
Due to the use of tensor-product basis and quadrature, the DG operators possess a Kronecker-product structure . Since the nodal points and quadrature points are collocated, there is no need for an interpolation operator. The one-dimensional mass matrix on the reference element, denoted , is given by a diagonal matrix with quadrature weights on the diagonal. The weighted one-dimensional differentiation matrix is obtained by evaluating the derivatives of the basis functions at the nodal points and multiplying by quadrature weights,
We will often make use of the following two simple properties of the differentiation matrix:
The weighted differentiation matrix satisfies the following two useful properties, which we will make use of extensively in this work:
On the reference element, the local mass and differentiation operators can be obtained through Kronecker products. For instance, for we have
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,
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 Gauss-Lobatto nodal points. These terms are evaluated using the freestream-preserving procedure descried in . 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 Gauss-Lobatto 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.
First, we consider conservation. The analogous statement at the discrete level is that
where is a vector of all ones. First, note that since the numerical flux function is single-valued, we have
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
proving the conservation property for the discretized weak form.
For the discretized strong form, we make use of property (P2). For any
and we therefore obtain conservation for the strong form as well:
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 left-hand side of (15) is zero. Therefore, the discretization will preserve constants if the following identity holds:
This so-called metric identity is discussed in detail in , and can be enforced through proper evaluation of the entries of . In this case, we have the following result:
Let satisfy the discretized strong form (15). Furthermore, suppose that is spatially constant. Then, .
3. Construction of the IDP low-order method
We now modify the discretization (15) with the goal of obtaining a method which is invariant domain preserving (IDP). In Guermond and Popov  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 high-order methods with large stencils results in very dissipative methods and typically poor-quality results, as observed in , 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 first-order accurate, the underlying sparse discretization is not required to be high-order accurate.
We make the following simple modification to the DG-SEM method described above. Notice that the wide stencil of the high-order method is a consequence of the fact that the one-dimensional differentiation matrix is dense. We therefore replace with a sparser version that is first-order accurate. is obtained by integrating the derivatives of piecewise linear basis functions on the mesh defined by the Gauss-Lobatto points in the interval . is therefore given by
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 column-sum properties (P1) and (P2). The modified operators then replace the standard DG-SEM derivative operators in the formulation (14), in order to obtain the following modified discretization:
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 zero-sum 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 high-order metric terms , but are designed such that the first-order discretization satisfies the metric identities. Since the modified low-order method is itself only first-order 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
where denotes the th component of the outward-facing normal at the Gauss-Lobatto points of face . This is equivalent to the following underdetermined system of equations:
Solving for the correction matrices , we obtain the system
Note that the differentiation matrix on the left-hand side has a null-space consisting of all constant functions. In order for a solution to exist, the right-hand side must be orthogonal to this null-space. In other words, the entries of the right-hand side vector must sum to zero. This can be seen to be true for the second term on the right-hand side by a simple application of the zero-sum property (P1). In order to show this property for the first term on the right-hand side, we make use of the summation-by-parts property (19):
) 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 right-hand 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 low-order discretization
As a consequence of the choice of metric terms, this discretization satisfies the following properties.
Let satisfy (29). Then, the following two properties hold:
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 Lax-Friedrichs 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 Lax-Friedrichs flux (equation (6)), the term of the form can be written as
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:
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
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
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 Lax-Friedrichs term may be combined with the graph viscosity by appropriately defining the viscosity coefficients. Let for be given by
The discretization given by (35) satisfies the following properties.
Since the low-order 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) Runge-Kutta method for the temporal discretization . 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
At this point we can introduce the so-called “bar states” (also referred to as “intermediate limiting states”, cf. ), which are defined by
Suppose is a convex invariant set of (1) such that Then, the bar states belong to .
See . ∎
Having defined the bar states, we rewrite (37) as a convex combination
As an immediate consequence of this rewriting is the following proposition.
Suppose the following CFL condition holds:
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
where is a prescribed velocity field. We assume that the velocity field is divergence free, i.e. . We are interested in ensuring that the low-order 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 low-order 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 . We define the modified bar states for the advection equation by
where , and denotes the coordinates of the th node of the mesh. The modified graph viscosity coefficients for the advection equation are given by
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.
The property was essential to writing the low-order 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 low-order IDP method to the graph viscosity method applied to the unsparisified DG-SEM 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 low-order IDP discretization (35) with the high-order discretization (15), to obtain a bound-preserving scheme. To this end, we consider a forward Euler time discretization of both the low-order and high-order approximations (as mentioned above, the extension to high-order time integration is straightforward using SSP methods). The provisional high-order update for element is given by
and the low-order IDP update is given by
We introduce the short-hand notation
Note that the provisional high-order solution may not satisfy convex invariants, such as maximum principles, positivity, etc. However, the low-order solution is guaranteed to be IDP as long as the CFL condition (40) is satisfied. We then define the flux-corrected update by
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 zero-sum property:
Let denote the set of indices associated with the element . Then,
In fact, the above proposition can be extended to consider “lines” of degrees of freedom in a dimension-by-dimension 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:
Consider a “line” of nodes within an element (for any and multi-index as described above). Then
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 low-order method is IDP, we have that if for all , then . We wish to enforce a similar property for the high-order 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.
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 low-order discretization:
4.1.2. Elementwise limiting for linear constraints
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 flux-corrected 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 one-dimensional 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 one-dimensional subcell limiting procedure along lines of nodes, in a dimension-by-dimension 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 “subcell-centered 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 :
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:
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,