An Entropy Stable High-Order Discontinuous Galerkin Method on Cut Meshes

While low-order methods can be extremely robust, their accuracy and efficiency per degree of freedom (DOF) can leave much to be desired. Low-order methods have higher numerical dispersion error than high-order equivalents [19] and are less efficient per degree of freedom [20]. As a result, high-order methods are increasingly preferred for fluid applications that are especially sensitive to diffusion, such as unsteady flows featuring vortices and turbulence [21], and the simulation over long distances or times [20]. However, the low-diffusion error of high-order methods is a double-edged sword: high-order methods tend to be unstable in the presence of shocks and discontinuities [20].

In the context of hyperbolic conservation laws, instability can be attributed in part to a method’s failure to satisfy the entropy inequality. Entropy stability can be imposed on a scheme externally by adding routines/operators that dissipate entropy to an otherwise entropy unstable scheme. Common approaches for externally enforcing entropy stability include artificial viscosity and regularization [22, 23], and filtering and/or limiting procedures [24, 25]. However, these approaches can result in the loss of high-order accuracy [20] and the need to tune parameters can result in a loss of robustness from one simulation to the next: too little dissipation and the scheme is unstable, too much and the solution is dominated by dissipation error.

A major difficulty in designing an inherently entropy stable high-order method is the use of the chain rule in space (as shown in (6)) in the derivation of the entropy statement as it has no semi-discrete equivalent. Summation-by-parts (SBP) operators paired with flux-differencing provides a means to side-step the need for the chain rule in the derivation of the semi-discrete entropy statement, a proof of which is presented in [26]. SBP operators are so named as they are a semi-discrete analog of integration-by-parts and as such allow information in a volume to be replaced by information on the volume’s boundary.

SBP operators have been established for a number of finite volume [10, 11, 14, 15, 12, 13] and DG schemes [27, 28, 29, 30]. SBP operators can also be constructed for finite difference schemes [31, 6]. SBP schemes are high-order accurate and can be used to create inherently entropy conservative/stable schemes.

Here, we use an entropy stable high-order SBP DG scheme, which pairs SBP operators with a flux-differencing DG formulation and entropy conservative/stable fluxes. In the DG setting, on a given element of the mesh classical diagonal-norm SBP schemes require the mass matrix $\bm{M}$ to be diagonal and positive-definite and SBP operators $\bm{Q}_{d}$ (which act as differentiation operators on volume of the element) and boundary integration matrices $\bm{B}_{d}$ satisfy the SBP property:

While entropy stable high-order SBP methods are well-established, their treatment is typically limited to simplicial and tensor product elements; in 2D spatial domains these correspond to triangular and quadrilateral elements, respectively. Here, we are interested in domains with embedded boundaries. Embedded boundaries are most commonly handled via unstructured, fitted triangular meshes but such meshes can require either $h$ adaptivity or curvilinear elements to resolve curved boundaries well without loss of high-order accuracy.

As high-order elements have more DOF per element than lower order, to maintain reasonable computational cost they are best used on relatively coarse meshes, where “relatively” is interpreted with respect to an adequately refined mesh for low-order versions of that method. Fitted triangular meshes can require a large amount of elements to resolve a boundary well, potentially leading to bloated computational costs regardless of order. Here, we instead use Cartesian cut meshes, which feature a simple, unfitted Cartesian background mesh “cut” by the embedded boundaries to yield a hybrid mesh of Cartesian and cut elements as shown in Figure 2.

Cut meshes have a long history in fluid dynamics problems, with early work in the 1970’s [32]. Cut meshes allow embedded boundaries to be represented with high accuracy while simultaneously maintaining a coarse mesh away from embedded boundaries. The cost of flux-differencing methods is driven by the number of evaluations of the numerical fluxes, which scales linearly with the number of faces in the mesh and nonlinearly with the polynomial degree of the solution. Cut meshes feature fewer faces and thus fewer evaluations of the numerical flux functions than an equivalent quad-tri mesh and (typically) fewer elements and faces than fitted triangular meshes, making them a computationally economic partner for high-order methods.

While cut meshes are attractive for use with high-order methods, they present a number of challenges. Namely, cut elements can be arbitrarily shaped and sized leading to potentially extremely restrictive CFL conditions (i.e., the small cell problem) and difficulty in constructing high-quality quadrature rules. High-quality quadrature rules are necessary for SBP operators, which classically require diagonal norm mass matrices among other conditions.

On non-hybrid meshes (i.e. those consisting of a single element type) the cost of constructing SBP-appropriate quadrature (which can be expensive to compute) and the corresponding SBP operators is mitigated by the ability to construct the operators once on a single reference element for all elements in the mesh. The arbitrary geometries of cut elements, besides requiring custom quadrature rules for every cut element, cannot be mapped to a single reference element and thus greatly increase the difficulty and computational cost of constructing traditional SBP operators on cut meshes.

Even when limited to known element geometries, SBP operators present a new issue when used on hybrid meshes: the conditions needed to induce the classical diagonal-norm SBP property can differ between different types of elements. For example, despite SBP operators being well-established on simplices and tensor-product elements individually, the differing conditions on each type of element create compatibility issues on hybrid quad-tri meshes [1]. However, relaxations of the SBP property exist, such as the generalized SBP property and hybridized SBP property [33, 26, 1], which reduce the cost of constructing SBP operators.

Of particular interest to us are skew-hybridized SBP operators, which were developed by Chan in [1] to overcome the compatibility issue on hybrid quad-tri meshes. Hybridized SBP operators were originally introduced as “decoupled SBP” operators in [26] but subsequently became known as “hybridized” SBP operators [34]. “Skew” hybridized SBP operators are so named due to their use of the skew-symmetric matrix $(\bm{Q}_{d}-\bm{Q}_{d}^{T})$. Given non-SBP differentiation operators $\bm{Q}_{d}$ and boundary integration operators $\bm{B}_{d}$, skew-hybridized SBP operators $\bm{Q}_{H,d}$ and $\bm{B}_{H,d}$ can be constructed via

and

for $d=1,2$. The skew-hybridized SBP operators satisfy the hybridized SBP property:

which is a block version of the traditional SBP property.

The matrix $\bm{E}$ used in (20) is an extrapolation operator taking the solution values at volume quadrature points and mapping them to values at surface quadrature points (its exact definition will be given later in Section 2.1.2). An important property of $\bm{E}$ is that it is dependent on both the volume and surface quadrature rules [26, 1]. With this dependence on volume and surface quadrature, skew-hybridized SBP operators can be thought of as a volume-based derivative operator ($\bm{Q}_{d}$) with boundary-based correction terms coming from the terms with $\bm{E}$ [1, 35].

Importantly, the accuracy of a skew-hybridized SBP operator as a differentiation operator is now tied to the accuracy of the both the surface and volume quadrature [1, 26]. However, since skew-HSBP operators do not require $\bm{Q}_{d}$ to be a SBP operator, they can be readily extended to cut meshes provided sufficiently exact quadrature rules can be constructed on cut elements.

A number of methods exist for constructing quadrature rules on cut meshes/arbitrary geometries. In some cases, the shape of cut elements is restricted (e.g. to simplicies), allowing established quadrature rules to be applied [36]. Approaches for constructing quadrature rules typically fall into two categories: generative methods, which numerically create a custom quadrature rule for a specific cut element, and subdivision approaches, which divide the cut element into a collection of subelements who have established quadrature rules (e.g., subtriangulation) from which a composite rule is constructed.

A common starting point for many generative approaches is moment-fitting, where quadrature rules are constructed on a set of points from target integrals [37, 38, 39]. Not all moment-fitting approaches guarantee polynomial exactness (which is often dependent on boundary representation) and positive quadrature weights, though some, such as [38, 40], do but may rely on using a possibly-extreme number of points to do so [41]. Such methods are also subject to the quality of the sampling points used as they affect the conditioning of the system.

There are however, other generative methods for creating quadrature rules. Work by Saye in [42, 43] provides an algorithm for constructing polynomially exact, positive weight quadrature rules for geometries defined by polynomial level-sets and is commonly used in the cut mesh community. However, Saye’s algorithm depends on geometries defined by polynomial level-sets, though the algorithm has also been extended to non-polynoimally defined level-sets via localized projection [44]. Here, however, we use explicit parameterizations to represent cut boundaries as we did previously in [45] and thus are not eligible to use Saye’s routines.

Subdivision approaches use existing quadrature on known geometries, such as simplicies and tensor-product elements, to construct a composite quadrature rule on a single cut element. Examples include quad/oct-tree quadrature [46, 36] and subtriangulation approaches [47, 48]. While subdivision approaches have the advantage of leveraging known quadrature formulas and efficient, well-developed libraries, they–by construction–contain many more quadrature points than is needed for a given polynomial degree by virtue of the subdivision. In the case of moving boundaries, the modularity of subdivision-based rules can be beneficial, however, for static boundaries the cost is less permissible.

While the skew-hybridized SBP formulation only requires polynomially exact, positive weight quadrature, the number of quadrature points is an important factor in the computational costs of a scheme both in runtime and memory. Therefore rules using an excessive number of points, such as non-negative least-squares and subdivision, are less preferable with respect to computational cost. However, schemes yielding/using a large number of points can also be quite robust and utilize well-establish libraries and thus warrant further consideration.

We construct our cut meshes using the Julia libraries PathIntersection.jl and StartUpDG.jl. Our cut meshes are constructed using explicit parameterizations of the embedded boundaries as previously detailed in [45]. Using explicit parameterization has a number of benefits, such as allowing piecewise-defined curves and thus an easy interface for imposing boundary conditions on portions of the cut boundaries and exact representation of sharp features. In particular, our code uses stop points to define both user-defined subintervals and junctions in piecewise curves.

Our code provides a piecewise representation of the boundary of every cut element, with stop points denoting the endpoints of the element’s faces. We use the stop points to triangulate each cut element via the Julia package Triangulate.jl. For a given cut element $D^{k}$, we decompose it into $N_{T}$ polynomially-defined non-overlapping curvilinear triangles $T_{i}$:

For each triangle, we construct a geometric mapping $\bm{g}_{d}:\hat{T}\to T_{i},\leavevmode\nobreak\ d=1,2$ from a reference triangle, $\hat{T}$, to each subtriangle $T_{i}$. Note for high-order accuracy the polynomial degree of the geometric mappings must be isogeometric, i.e., at least the same degree as the degree of the solution. We take the mappings to be total degree $N$ in all meshes, i.e. $\bm{g}_{d}\in\mathbb{P}^{N}(\hat{T})$. An example of the subtriangulation process is shown in Figure 3.

The geometric mappings provide both a polynomial approximation of the typically non-linear curved boundaries and access to a reference element, allowing both exact surface and volume quadrature rules to be constructed. However, the geometric mapping and its Jacobian also increase the degree of the integrand.

We will illustrate the effect of the geometric mapping with respect to the volume integrals. Consider a polynomial $p\in\mathbb{P}^{M}(D^{k})$. On a given subtriangle $T_{i}\subseteq D^{k}$, to map an integral with respect to $\bm{x}\in T_{i}$ to an integral with respect to $\hat{\bm{x}}\in\hat{T}$, we take $\bm{x}_{d}=\bm{g}_{d}(\hat{\bm{x}})$, $d=1,2$. The subintegral on $T_{i}$ is then related to its corresponding integral on $\hat{T}$ by

For $p\in\mathbb{P}^{M}(D^{k})$ and $\bm{g}_{d}\in\mathbb{P}^{N}(\hat{T}),\leavevmode\nobreak\ d=1,2$ we have that $p\circ\bm{g}\in\mathbb{P}^{MN}(\hat{T})$. Similarly, the determinant of the Jacobian $\left|\nicefrac{{\partial\bm{g}}}{{\partial\hat{\bm{x}}}}\right|\in\mathbb{P}^{2(N-1)}$. Thus the final integrand on the reference triangle is a polynomial of total degree $MN+2(N-1)$. For the accuracy conditions on $\bm{Q}_{H,d}$, we need $M=2N-1$, resulting in the final integrand $\hat{p}$ on the reference triangle being total degree $(2N^{2}+N-2)$; a quadratic increase in degree with respect to $N$, the polynomial degree of the solution. To exactly integrate the mass matrix would result in the final integrand being degree $2(N^{2}+N-1)$. Table 1 shows the polynomial degree of the final integrand on the reference triangle for various values of $N$ for reference.

Table 1 helps to illustrate why subdivision type quadrature approaches–when used by themselves–are so inefficient when paired with high-order methods: the subquadrature rules must be constructed for polynomial degrees far in excess of the degree of the solution. Moreover, the cost of these many point subquadrature rules is multiplied by the number of triangles in the triangulation of each cut element.

However, subdivision based quadrature rules are good starting point for a pruning routine: by construction the subquadrature rules are exact and have purely positive weights. We use the QR factorization approach of van den Bos [51] for pruning, which we describe in Appendix 6.

While not the focus of this paper, one of the greatest challenges of/deterrents from using cut meshes is the small cell problem: cut elements can be arbitrarily small/skewed. When using explicit time integration schemes, the small cell problem can result in extremely prohibitive CFL conditions. One method for addressing the small cell problem is state redistribution, which was introduced [58] and adapted to high-order methods in [59].

State redistribution is high-order accurate and conserves the average of the solution. State redistribution involves projections spanning multiple elements and averaging of those projections (or more generally convex combinations of those projections [60, 61]) and as a result is dissipative in nature. In the context of linear symmetric hyperbolic conservation laws, state redistribution can be added to a $L_{2}$ energy stable scheme without damaging the scheme’s energy stability/conservation as shown in [45]. However, state redistribution readily violates conservation of entropy when applied to an entropy conservative scheme.

Encouragingly, in the entropy stable case we have observed that state redistribution does not violate the entropy stability of our scheme. We show the effect of state redistribution on a system’s entropy in Section 3.2. In keeping with our observations, many of our entropy stable numerical experiment use state redistribution. However, we acknowledge that it remains to be proven whether state redistribution can in fact violate entropy stability; we save this investigation for future work.

In the following numerical experiments we consider two hyperbolic systems, the shallow water equations and the compressible Euler equations, on 2D domains. The definitions of these systems and their respective fluxes and the time integration scheme used are given in this section for reference.

In both systems, when an entropy stable flux is needed we use the Lax-Friedrichs flux

where $\lambda$ is the maximum wavespeed or an estimate thereof as previously mentioned. We use the Lax-Friedrichs fluxes provided in the Trixi.jl library, which use the wavespeed approximation by Davis [64].

Here we numerically verify the entropy status of our scheme. For this experiment, we are interested in the evolution of the system’s entropy over time. This can be done by considering the integral over the domain of the time derivative of entropy which, for entropy conservative boundary conditions, satisfies

with equality holding in the absence of shocks. At the semi-discrete level, this entity is called the entropy residual and is given by:

where the superscript $(k)$ indicates the appropriate entity on the $k^{th}$ element and

In the semi-discrete case, equality (up to round-off error) holds if the scheme is entropy conservative, while the inequality hold for an entropy stable scheme.

We first consider the shallow water equations on a Cartesian embedding domain $\hat{\Omega}=[-1,1]^{2}$ with a circle of radius $R=0.331$ removed from the center. We impose a $16\times 16$ background Cartesian mesh. To test the entropy status of our scheme we compute the entropy residual as shown in (56) at each time step. We use an initial time step of $\Delta t=1\times 10^{-4}$ and integrate up to an end time of $t=1$.

We enforce reflective wall boundary conditions, which are entropy conservative [69], on all boundaries including the circle walls. We start with a zero-velocity solution with discontinuous water height given by

Note that in the continuous setting, weak solutions should dissipate entropy thanks to the discontinuity in the water height. This serves as a stress test for an entropy conservative scheme: it should conserve entropy despite the discontinuity, which in non-entropy conservative/stable scheme can lead to instability. In the entropy stable case, entropy should immediately be dissipated due to the discontinuity and the dissipation wane as the solution approaches its steady state. Figure 4 shows the initial water height.

We consider four schemes for this experiment: the entropy conservative scheme with and without state redistribution, and the entropy stable scheme with and without state redistribution. All schemes use polynomial degree $N=4$. Figures 5 and 6 show the entropy residual for each of the four schemes.

Without state redistribution, the entropy conservative scheme maintains entropy conservation within round-off error despite the presence of the shock. The addition of state redistribution to the entropy conservative scheme however violates entropy conservation by a non-negligible margin. The entropy stable scheme meanwhile maintains entropy stability even when state redistribution is applied. While it is encouraging to see that state redistribution did not violate entropy stability in this instance, we again stress that entropy stability may not always hold under state redistribution.

For this experiment we numerically verify the high-order accuracy of our scheme via an $h$-convergence study. We take the Cartesian embedding domain $\hat{\Omega}=[-1,1]^{2}$ from which we remove a circle of radius $R=0.331$ from the origin. We consider the manufactured solution to the shallow water equations defined by primitive variables:

Figure 7 shows the water height at time 0. We start with a $4\times 4$ Cartesian background mesh and refine in each dimension by a factor of 2 up to a $32$ background mesh. We use the manufactured solution values as boundary conditions and simulate up to an end time of $t=0.3$. We use an initial time step of $\Delta t=1\times 10^{-4}$ for the $4\times 4$ mesh and $\Delta t=1\times 10^{-5}$ for the other meshes.

We conducted the convergence study for solutions of degree $N=2,3,4$. As shown in Figure 8, our method recovers the expected convergence rate of $h^{N+1}$ in each case.

Here we consider an analytic solution to the compressible Euler equations. When initialized with constant velocity and pressure, the compressible Euler equations reduce to the advection equation with vector-valued coefficients applied to the initial density field. This special instance of the Euler equations is known as an entropy wave and is one of three modes of linear flow fluctuation in gases (the others being acoustic waves and vorticity waves) [70]. The reduced equations are:

which can also be expressed as

For this experiment, we again consider the embedding domain $\hat{\Omega}=[-1,1]^{2}$ over which we impose a $16\times 16$ background Cartesian mesh. A circle of radius $R=0.331$ centered at the origin is removed. We consider the initial condition

which by the method of characteristics applied to the reduced PDE/advection equation (60) yields the analytic solution

For boundary conditions, we prescribe the analytic solution on all boundaries. For time integration we use an initial time step of $1\times 10^{-6}$ (a typical timestep was on the order of $1\times 10^{-4}$) with state redistribution. The simulation was run to an end time of $t=1.3$ with polynomial degree $N=4$. Figure 9 shows snapshots the solution an error at various points in time alongside the error in the density field. As illustrated in the error plots, cut elements produce higher errors than Cartesian elements, however, the maximum error magnitude is five orders of magnitude less than the magnitude of the solution.

For this experiment we consider the compressible Euler equations to simulate supersonic flow over a biconvex airfoil. This experiment highlights our method’s ability to handle sharp features in the mesh and solution while remaining stable.

For the mesh we take the Cartesian embedding domain $\hat{\Omega}=[-0.5,5.5]\times[-1.5,1.5]$ over which we impose a $120\times 60$ background Cartesian mesh. The airfoil, whose unscaled geometry is given in Figure 10, is scaled to a non-dimensionalized chord length of 0.5 and centered about the $x$-axis at zero angle of attack and with its chord line at $y=0$.

We apply reflective wall boundary conditions to the airfoil and the top ($y=1.5$) and bottom ($y=-1.5$) walls of the domain. We take the left ($x=-0.5$) boundary as the inlet with density $\rho_{\infty}=50$, pressure $p_{\infty}=50$, $x$-velocity $u_{\infty,1}=$ Mach 1.5, and $y$-velocity $u_{\infty,2}=0$. Extrapolation boundary conditions were used on the right ($x=5.5$) boundary.

For the initial condition we use an impulsive start with the solution initialized to the freestream conditions. We take $N=4$ and simulate the system up to time $t=5$ with state redistribution using an initial time step of $\Delta t=1\times 10^{-8}$ (a typical timestep was on the order of $1\times 10^{-4}$). During this time period, a bow shock immediately forms off the leading edge of the airfoil before reflecting off the top and bottom walls as shown in Figure 11. Our method is able to capture these shock waves and shock-shock interactions without loss of stability.

For our last experiment we consider the shallow water equations on a domain with three embedded objects. We take the Cartesian embedding domain to be $\hat{\Omega}=[-1,1]\times[0,5]$ over which we impose a $16\times 40$ background Cartesian mesh. Figure 12 shows the mesh for this problem. The three embedded circles share the same radius $R=0.31$ and are centered at the points $(x,y)=(0,1),(0,2.5),(0.4)$. We again use a zero velocity initial condition and take the initial water height to be the discontinuous function

We impose reflective wall boundary conditions on all boundaries and take the solution degree as $N=4$. We simulate the system up to an end time of $t=5$ using an initial time step of $\Delta t=1\times 10^{-4}$ with state redistribution. This end time is sufficient for the front to propagate to and reflect off of the far wall. Figure 13 shows snapshots of the solution at various points in time.