Improvements of algebraic flux-correction schemes based on Bernstein finite elements

2.1 Model problem

Let Ω ⊂ R d , d ∈ {1, 2, 3}, x ∈ Ω, t ⩾ 0, and let u ( x , t ) ∈ R M , M ∈ N , be the solution to the conservation law

where ∇⋅ is the matrix-divergence operator applied to the inviscid flux f = f ( u ) ∈ R M × d . To formulate a well-posed initial-boundary value problem, (1) has to be equipped with the initial condition u(⋅, 0) = u ₀ and suitable data u 0 = u 0 ( x ) ∈ R M as well as appropriate boundary conditions. The latter are specified for the test problems in Section 6. For many specific applications (e.g., scalar problems, shallow water equations, gas dynamics), the exact solution to (1) is known to satisfy certain admissibility constraints that can be described by constraining u to lie in a some convex set A ⊂ R M . For instance, in the scalar case, M = 1, A is an interval bounded by the extrema of initial and boundary conditions [36], Ch. 6]. Suppose there exists a convex function U : A → R and a corresponding flux F : A → R d such that F ′(u) = U′(u) ^⊤ f ′(u), then (U, F ) is called an entropy pair for the conservation law (1). By the chain rule, a conservation law for the entropy U can be derived under the assumption that u and U are smooth. In general however, weak solutions to the nonlinear problem (1) can develop discontinuities in finite time. Using, for instance, the concept of vanishing viscosity solutions [36], Sect. 6.3], one can show well-posedness of weak entropy solutions for the scalar case, if a weak form of the entropy inequality

holds for all (U, F ) that are entropy pairs for (1). If discrete versions of (2) are satisfied, ideally in a localized manner, the occurrence of nonphysical weak solutions can be averted in numerical approximations, which motivates the use of entropy-stable approximations, e.g., [3], [4], [7], [12], [21].

2.2 Generic monolithic convex limiting discretization

We now summarize the concept of monolithic convex limiting (MCL), first proposed by Kuzmin [13], see also [7], [18], [21], [25], [32]. This algebraic flux correction (AFC) paradigm operates on the semi-discrete level and is capable of enforcing a variety of constraints. A generic MCL semi-discretization can be written as

Here the index i ∈ {1, …, N} refers to a node, N is the number of unknowns for each variable, and thus the total number of unknowns is MN. Note that in the system case all components of the solution vector u are interpolated on the same set of nodes, which are associated with a certain node x i ∈ Ω ̄ , i.e., we do not consider staggered approaches. Similarly, m _i > 0 is associated with this node and typically represents an entry of a diagonalized mass matrix. The set N i ⊂ { 1 , … , N } is the nodal stencil of x _i , i.e., the set of nodes x _j that directly interact with x _i . In piecewise linear continuous finite element methods, N i is the set of mesh vertices that are nearest neighbors to node x _i in the sense that they are endpoints of a certain mesh edge. Entities featuring two distinct indices such as the symmetric artificial diffusion coefficient d _ij = d _ji ⩾ 0 refer to terms depending on the two nodal values u _i and u _j . In the MCL framework, the flux-corrected bar states u ̄ i j * [13] are algebraically enforced to satisfy user-defined constraints of convex nature. Upon temporal discretization of (3) with a strong-stability preserving (SSP) Runge–Kutta (RK) method [37], [38], [39], Shu–Osher updates can be written as a convex combination of forward Euler steps

where τ > 0 is the time step. Suppose that u i ∈ A i for a certain convex subset A i ⊂ A of the largest admissible set. If u ̄ i j * ∈ A i can be guaranteed for j ∈ N i \ { i } and the CFL-like time step restriction

holds, then u i F E is a convex combination of the solution at the previous iteration u _i and the flux-corrected bar states u ̄ i j * . In conclusion, by convexity of A i , we then have u i F E ∈ A i as well. Recent works, e.g., [33], [40], also demonstrate how AFC and MCL can be performed using non-SSP time discretizations, including implicit methods and applications to viscous problems. Although they appear promising, these issues shall not be addressed here since our focus lies on the spatial discretization. The particular choice of the sets A i and how to enforce u ̄ i j * ∈ A i for all j ∈ N i \ { i } sets different MCL techniques apart from each other. Before discussing high-order DG discretizations, we give a brief example of arguably the simplest MCL variant.

2.3 Case study: MCL based on classical finite elements for the Euler equations

We only summarize the steady case here. More information can be found in ref. [7] and the references therein. Let I be the d × d-identity matrix and let ρ = ρ ( x , t ) ∈ R , v = v ( x , t ) ∈ R d , p = p ( x , t ) ∈ R denote density, velocity, and pressure of an ideal polytropic gas. We consider the compressible Euler equations

where u = [ ρ , ρ v ⊤ , ρ E ] ⊤ is the vector of conserved unknowns (density, momentum, and total energy), and

where γ > 1 is the adiabatic constant and e is the specific internal energy. The (physically motivated) largest admissible set for (6) and (7) is A = { u = [ ρ , ρ v ⊤ , ρ E ] ⊤ ∈ R d + 2 : ρ ⩾ 0 , p ⩾ 0 } . Note that due to (7), nonnegativity of pressure is equivalent to the physically motivated constraint that internal energy e may not become negative. The usual entropy pair for (6) and (7) is given by (U, F ) with U = −ρlog(pρ ^−γ ), F = v U [41].

Let us discretize (1) using a conforming simplicial mesh and local P 1 -polynomials with corresponding Lagrange basis functions φ i ∈ C ( Ω ̄ ) , i ∈ {1, …, N}, defined implicitly by the interpolation property φ _i ( x _j ) = δ _ij for all i, j ∈ {1, …, N}. We then set m _i = ∫_Ω φ _i  d x , c _ij = ∫_Ω φ _i ∇φ _j  d x and

For a unit vector n _ij = c _ij /| c _ij |, (8) corresponds to the largest absolute eigenvalue of the Jacobian of f (u _i ) or f (u _j ) projected onto n _ij . This estimate for the exact wave speed does not produce states outside of the set A [42], App.], which is why we use (8) to define the diffusion coefficients as follows [24], [43], [44]:

Next, we define the low-order bar states u ̄ i j = u i + u j 2 + ( f i − f j ) c i j 2 d i j [23], [24] and their MCL counterparts [13]:

where f _i := f (u _i ) and f i j * = − f j i * is a flux-corrected counterpart of the target flux f _ij , which for steady-state calculations can simply be set to f _ij = d _ij (u _i − u _j ). The low-order counterpart u ̄ i j of u ̄ i j * initially appeared in Harten et al. [23]. Guermond and Popov [24] were the first to use it in the AFC context to analyze the low-order method that is the first-order predictor in FCT-type convex limiting schemes and, similarly, the low-order version of MCL approaches. In either case, the limited antidiffusive fluxes f i j * are set to zero.

In general, the least-restrictive constraint that should be enforced for any hyperbolic problem is the admissibility of solutions, i.e., u ∈ A . For the compressible Euler equations, we need to ensure that densities remain nonnegative, which can be achieved by enforcing positivity of the first component of u ̄ i j * = ρ ̄ i j * , ( ρ v ) ̄ i j * , ( ρ E ) ̄ i j * . Using definition (10) and skew symmetry of f i j * , this task is accomplished by setting the first component of f i j * to max − 2 d i j ρ ̄ i j , min { f i j ρ , 2 d i j ρ ̄ j i } , where ρ ̄ i j and f i j ρ correspond to the first component of the low-order bar state u ̄ i j and the unlimited flux f _ij , respectively. Next, we can enforce additional constraints such as nonnegativity of pressure and Tadmor’s entropy stability condition [30], [31] by further reducing the magnitude of f i j * . These issues are described in detail in refs. [13], [18], [27], respectively, see also [7]. Additionally, the quality of numerical solutions can be improved by enforcing discrete maximum principles in addition to nonnegativity, see, e.g., [7], [10], [16], [21], [28], [40]. All of these tasks (and more, for instance, well-balancedness in the case of the shallow water equations [12]) can be guaranteed simultaneously.

Since nonnegativity of pressure is a nonnegotiable constraint (as is nonnegativity of density), we briefly summarize the typical MCL pressure fix [13], [21], [25]. For a pair of nodes i ≠ j let f i j * be the prelimited antidiffusive flux. Setting f i j * * = α i j f i j * , where the correction factor α _ij = α _ji ∈ [0, 1] is defined via [13]:

where w ̄ i j = 2 d i j u ̄ i j = w ̄ i j ρ , w ̄ i j ρ v ⊤ , w ̄ i j ρ E ⊤ , f i j * = f i j ρ , * , f i j ρ v , * ⊤ , f i j ρ E , * ⊤ , and

guarantees that the pressure of the bar state u ̄ i j * * = u ̄ i j + f i j * * / ( 2 d i j ) remains nonnegative. The derivation of (11) can be found in refs. [13], [28]. The following result is an interesting observation made in the process of this work.

3.1 General considerations and second-order continuous Galerkin schemes

Let us once more consider a P 1 or Q 1 continuous Galerkin discretization as described in Section 2.3. If the domain has only periodic boundaries, then each component C ^k of the discrete gradient operator C = ( C 1 , … , C d ) = ( c i j ) i , j = 1 N is a skew-symmetric N × N matrix due to the divergence theorem. Otherwise,

where n = n ( x ) ∈ R d denotes the unit outward normal to ∂Ω. In this setting, each matrix C ^k , k ∈ {1, …, d}, is almost skew-symmetric, with only a few pairs of entries where b _ij ≠ 0 being exceptions to this rule. For flux-correction schemes, the skew symmetry of C is often desirable to design new algorithms [33], simplify computations (note that c _ij = − c _ji implies λ _ij = λ _ji in (8) (the computation of λ _ij is a significant cost factor in AFC schemes for systems), and to prove certain theoretical results such as entropy stability [18], see also Section 4.3.2. Therefore, it is unfortunate that a few boundary nodes invalidate skew symmetry.

We now explain how the issue raised at the beginning of this section can be resolved starting in the already introduced CG setting. Let us take a step back and remind ourselves how the matrices C arise in the first place. The strong form of the CG discretization of the conservation law (1) contains the integral

which can generally only be approximated via quadrature. The group finite element formulation [45], [46]:

interpolates the inviscid flux f (u _h ) in the finite element space. On simplicial meshes with periodic boundary conditions, this approach can be interpreted as nodal quadrature [46]. Inserting (13) into (12), we obtain

which is how the discrete gradient operators arise in AFC schemes. The last equality in (14) holds because the row sums of each C ^k are zero due to the partition-of-unity property ∑ j = 1 N φ j ( x ) = 1 for any  x ∈ Ω ̄ . Replacing c _ij in (14) by some c ̃ i j ∈ R d such that c ̃ i j = 0 if j ∉ N i is equivalent to adding

to (14). If for all k ∈ {1, …, d} the row sums of each C ̃ k in C ̃ = ( C ̃ 1 , … , C ̃ d ) = ( c ̃ i j ) i , j = 1 N sum to zero, and the columns of the matrices C ^k and C ̃ k add up to the same values c j k stored in c j = c j 1 , … , c j d , then

Thus, the proposed modification does not lead to a change in the global mass balance. Note that by writing the correction term (15) using a sum over only indices j ≠ i, we may choose the diagonal entries c ̃ i i arbitrarily. This possibility is what will allow us to define the skew-symmetric discrete gradients below.

An additional constraint that is desirable in the context of high-order AFC schemes is the sparsity of C ̃ [10], [14], [20], [21], [25], [32], [47]. By that we refer to the property that every nonzero offdiagonal entry corresponds to two distinct nodes that are closest neighbors (in a certain metric) in the submesh obtained by using each Bernstein node as a degree of freedom of a continuous subcell P 1 or Q 1 (depending on the element geometry) discretization [13], [14], [25]. The corresponding stencil N ̃ i will be further specified below. Each pair of nodes (i, j) with j ∈ N ̃ i \ { i } corresponds to one antidiffusive flux that needs to be limited and incorporated. Thus, the fewer elements in the stencils N ̃ i the better in terms of computational resources.

Thus, we may replace the discrete gradient operator C by matrices C ̃ of the same size satisfying

The remainder of this section is dedicated to constructing such matrices C ̃ . For second-order continuous Galerkin schemes with weakly-enforced boundary conditions, C satisfies (16b)–(16d). Thus finding a C ̃ that also fulfills (16a) is simple (which is why no originality is claimed here, although we found no works using these matrices in a similar context). Given the original operators C and symmetric boundary matrices B = ( B 1 , … , B d ) = ( b i j ) i , j = 1 N , we only have to adjust a few entries of C , which can be done as follows:

Proving the validity of (16) for this operator is straightforward. So far, we have focused on global matrices. Below, we discuss elementwise operators because these become relevant in the high-order and/or DG cases.

3.2 Skew-symmetric discrete gradient operators for Bernstein polynomials

In Section 3.1, we established that the group finite element formulation [45], [46] produces the discrete gradients

where { φ i } i = 1 N are the Lagrange basis functions of P 1 or Q 1 continuous finite element spaces. Let us now consider elementwise operators C for a certain reference element K ⊂ R d and the Bernstein basis. These matrices are obtained by replacing the number of rows and columns N and the integration domain Ω in (17) by the local number of unknowns and the reference element K, respectively. With some abuse of notation, we avoid making this distinction, to avoid redefining all operators or carry additional indices. For each of the below element types, we construct matrices C ̃ = ( C ̃ 1 , … , C ̃ d ) = ( c ̃ i j ) i , j = 1 N , satisfying conditions (16).

3.2.3 Simplices

We now move on to the reference simplex K = conv { 0 d , e 1 , … , e d } ⊂ R d , where conv{} denotes the convex hull and e _i is the ith unit vector in R d . The barycentric coordinates are Λ 0 ( x ) = 1 − ∑ i = 1 d x i , Λ _i ( x ) = x _i for i ∈ {1, …, d}, and the Bernstein nodes shall be numbered using a multi-index notation i = (i ₀, …, i _d ) similar to the case of box elements. The set of multi-indices representing nodes in K is

J = J ( p , d ) = ( i 0 , … , i d ) ∈ { 0 , … , p } d + 1 : ∑ j = 0 d i j = N : = p + d p

and the (isotropic) pth-degree Bernstein basis functions read φ k p ( x ) = p ! ∏ l = 0 d k l ! ∏ l = 0 d Λ l ( x ) k l , k ∈ J .

3.2.3.1 The P 1 triangle

Before studying the general simplex case, let us briefly focus on the P 1 -triangle in 2D. By (16a)–(16c), all matrices satisfying (16) can be expressed as

C 1 = 1 4 − 1 a 1 − a − a 1 a − 1 a − 1 1 − a 0 , C 2 = 1 4 − 1 b 1 − b − b 0 b b − 1 − b 1 , a , b ∈ R .

The magnitude of nonzero off-diagonal entries in these matrices affects the time step restriction of convex limiting schemes due to (9) and (5). Thus, all nonzero off-diagonal entries of C ¹ and C ² should generally have the same absolute value if possible (see also Section 5 for further discussions on this issue). By this argument, it makes sense to only choose a, b ∈ {0, 0.5, 1}. Figure 1 shows various choices including three (Figure 1b–d) out of the possible nine combinations of parameters a and b in addition to the usual, non-skew-symmetric matrices C (Figure 1a) along with a variant, where a, b ∉ {0, 0.5, 1} (Figure 1e).

Figure 1:

Five different options for how to choose C ̃ on triangles, red represents entries in C ¹, blue represents entries in C ².

The case a = 1, b = 0 leads to dimensional decoupling and inhibits optimal convergence rates, see the numerical examples in Section 6.2. Setting a = 0, b = 1 corresponds to an unusual gradient, where coupling occurs in a direction that is opposite to the one in the previous case. This scheme exhibits optimal convergence behavior but requires significantly more time steps than the choice a = b = 0.5, which seems to be optimal in terms of CFL conditions. The gradient in Figure 1e corresponds to a discrete gradient that satisfies conditions (16) but uses off-diagonal entries of varying magnitude. As a result, the resulting CFL condition is much more restrictive, however, this method also exhibits optimal convergence rates. We have not experimented with gradients other than the ones shown in Figure 1 because the two matrices C ¹, C ² corresponding to any other choice of a, b ∈ {0, 0.5, 1} would be dimensionally inconsistent to each other.

In conclusion, uniqueness of simplicial discrete gradients satisfying conditions (16) does not hold. Note that for box elements the operators specified earlier are only unique because of their tensor-product structure (not possible here). Based on the summarized results obtained with each of the schemes in Figure 1 (detailed in Section 6.2), it makes sense to generalize the gradient in Figure 1d. In the remainder of this section, we demonstrate that this approach is feasible for general high-order simplices in arbitrary space dimensions.

3.2.3.2 General simplices

The realization that the gradient in Figure 1d appears to be optimal for P 1 -triangles implies the following structure for arbitrary p , d ∈ N : Conditions (16a)–(16c) fully determine the values of diagonal entries since

c j k = ∑ i = 1 N c ̃ i j k = ∑ i = 1 N c ̃ i j k + ∑ i = 1 N c ̃ j i k = ∑ i = 1 N c ̃ i j k + c ̃ j i k = 2 c ̃ j j k + ∑ i = 1 i ≠ j N ( c ̃ i j k − c ̃ i j k ) = 2 c ̃ j j k ⇒ c ̃ j j k = c j k 2 .

In any matrix row i, the offdiagonal entry in the jth column is nonzero only if the multiindices i and j correspond to nodes that are nearest neighbors. The magnitude of offdiagonal entries should be constant, and their sign should be determined by the spatial direction. If the node x _j lies in positive direction w.r.t. any Cartesian coordinate, then c ̃ i j k > 0 and c ̃ i j k < 0 otherwise for all k ∈ {1, …, d}, i.e., for all matrices. For closest neighbors other than these nodes (diagonal neighbors), the sign of c ̃ i j k depends on whether the diagonal direction is positive w.r.t. the coordinate x _k . If d > 2, there are diagonal neighbors with the same multiindex component for the x _k direction. The corresponding indices are set to zero in C ̃ k .

To formalize these considerations mathematically, we require some notation. For consistency with previous works [25], [32], we define connectivity on simplices as follows.

Definition 1.

Let i , j ∈ J . If ∃k, l ∈ {0, …, p}, k ≠ l, such that j = i + e _k − e _l , where e _k,l are the kth and lth unit vectors in R d + 1 , respectively, then we say j ∈ N ̃ i , i.e., j is in the local stencil N ̃ i .

In addition, let

(21) s = 1 2 ( d − 1 ) ! p + d − 1 p = p ! 2 ( p + d − 1 ) !

denote half the mass of a pth-degree Bernstein polynomial on the (d − 1)-dimensional reference simplex. Furthermore, let c j k = ∫ ∂ K φ j n k d s be as before, and for i , j ∈ J , define

(22) c ̃ i j k : = − s / d     if  ( j ∈ N ̃ i ) ∧ [ ( j 0 = i 0 + 1 ) ∨ ( j k = i k − 1 ) ] , s / d     if  ( j ∈ N ̃ i ) ∧ [ ( j 0 = i 0 − 1 ) ∨ ( j k = i k + 1 ) ] , − s     if  ( i = j ) ∧ ( i 0 > 0 ) ∧ ( i k = 0 ) , s     if  ( i = j ) ∧ ( i 0 = 0 ) ∧ ( i k > 0 ) , 0     otherwise.

Lemma 5.

The matrices C ̃ k = ( c ̃ i j k ) i , j ∈ J = ( c ̃ ( i 0 , … , i d ) ( j 0 , … , j d ) k ) i , j ∈ J given by (22) satisfy conditions (16).

Proof.

By Definition 1, i ∉ N ̃ i for all i ∈ J . Therefore, all cases in (22) are exclusive, and c ̃ i j k is well defined. Sparsity (16d) is built into Definition 1, and skew symmetry (16a) is proved by observing that for i ≠ j,

c ̃ j i k = − s / d     if  ( i ∈ N ̃ j ) ∧ [ ( i 0 = j 0 + 1 ) ∨ ( i k = j k − 1 ) ] , s / d     if  ( i ∈ N ̃ j ) ∧ [ ( i 0 = j 0 − 1 ) ∨ ( i k = j k + 1 ) ] , 0     otherwise, = s / d   if  ( j ∈ N ̃ i ) ∧ [ ( j 0 = i 0 + 1 ) ∨ ( j k = i k − 1 ) ] − s / d   if  ( j ∈ N ̃ i ) ∧ [ ( j 0 = i 0 − 1 ) ∨ ( j k = i k + 1 ) ] 0   otherwise = − c ̃ i j k .

We show (16c) by proving that

(23) c ̃ j j k = 1 2 c j k ∀ j ∈ { 1 , … , N } , k ∈ { 1 , … , d } ,

which, together with (16a) and (16b) (to be proven last), implies (16c). To better illustrate how (23) can be shown, let us consider the example sketched in Figure 2. The red nodes are precisely those corresponding to the nonzero diagonal entries in (22). These are set in accordance with (23), as are the entries for black and blue nodes: Diagonal entries for the former are zero because these nodes are either within the element interior or lie on a boundary where the respective component of the normal appearing in c j k is zero. The blue nodes are degrees of freedom for which the Bernstein polynomial is nonzero on more than one boundary segment, with the corresponding normal component being nonzero on precisely two (also for d > 2) of these, which depends on k. These two integrals cancel, which is consistent with definition (22).

Let us now rigorously prove (23) for (22) by formalizing these exemplified considerations in all required cases: If j ₀ j _k > 0 (black nodes), then c j k = 0 and (22) is consistent with (23) because c ̃ j j k = 0 . For j ₀ > 0 = j _k (red nodes not opposite the origin), the integral c j k = ∫ ∂ K φ j n k d s can be reduced to a (d − 1)-dimensional reference simplex, on which φ _j is a Bernstein polynomial in d − 1 variables. The component n _k of the normal to this boundary is always −1, thus, c j k = − 2 s by (21), and c ̃ j j k = − s . If j ₀ = 0 and j _k > 0 (red nodes opposite the origin), the situation is exactly reversed, and by similar arguments (including transformation to the reference simplex), we obtain c j k = 2 s = 2 c ̃ j j k . In the case j ₀ = j _k = 0 (blue nodes), contributions to the integral c j k arise from precisely two boundaries, one of which is always opposite the origin. These two integrals are clearly of opposite sign. Invoking transformation rules, one can show that their magnitudes are equal. Hence, c j k = 0 , which is consistent with (22) because c ̃ j j k = 0 . Thus, we have shown (23).

It remains to prove (16b), where, for clarity, we also distinguish between all relevant cases based on the row multiindex i ∈ J : For i ₀ = p, we have c ̃ i i k = − s , precisely d positive, and no negative entries in the row. Similarly, for i _k = p, c ̃ i i k = s , and the other d nonzero row entries are negative. For any other vertex, i.e., ∃l ∈ {1, …, d}\{k} with i _l = p (for d = 2, these are blue nodes in Figure 2), we have c ̃ i i k = 0 . The nodes with j ₀ = 1, j _l = p − 1 and j _k = 1, j _l = p − 1 each contribute one negative and one positive row entry, of which there are no more nonzero ones. We have now dealt with all corners of the simplex. At this stage, it makes sense to restrict ourselves to the case d > 1 (in 1D, it can easily be shown that (22) coincides with (18)). We define q := |{l ∈ {1, …, d}\{k} : i _l > 0}|, where |⋅| denotes the cardinality of a set. For i ₀ = i _k = 0, we have c ̃ i i k = 0 and there are precisely q negative and q positive entries in row i, which correspond to j ₀ = i ₀ + 1 and j _k = i _k + 1, respectively. For i ₀ = 0 < i _k < p, we have c ̃ i i k = s . The negative row entries correspond to either j = i + e₀ − e _k (one entry), j ₀ = i ₀ + 1 (q additional entries), or j _k = i _k − 1 (d − 1 additional (diagonal) neighbors). Moreover, there are precisely q positive entries in that row, which correspond to j _k = i _k + 1. Thus, the row sum is s − s d ( 1 + q + d − 1 ) + s d q = 0 . The case i _k = 0 < i ₀ < p is similar, and we obtain − s − s d q + s d ( 1 + d − 1 + q ) = 0 for the row sum. Finally, for 0 < i _0,k < p, we have c ̃ i i k = 0 and the row sum is − s d ( 1 + q + d − 1 ) + s d ( 1 + d − 1 + q ) = 0 . We have thus shown (16b), which together with (23) concludes the proof of (16c) as well as that of Lemma 5. □

Figure 2:

Bernstein nodes of the P 4 triangle and corresponding multiindices.