{'fraction_non_alphanumeric': 0.030566420396928872, 'fraction_numerical': 0.010243993294840752, 'mean_word_length': 4.923265506251532, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Positivity-Preserving Entropy-Based Adaptive Filtering for Discontinuous Spectral Element Methods

In this work, we present a positivity-preserving entropy-based adaptive filtering method for shock capturing in discontinuous spectral element methods. By adapting the filter strength to enforce positivity and a local discrete minimum entropy principle, the resulting approach can robustly resolve strong discontinuities with sub-element resolution, does not require problem-dependent parameter tuning, and can be easily implemented on general unstructured meshes with relatively low computational cost. The efficacy of the approach is shown in numerical experiments on hyperbolic and mixed hyperbolic-parabolic conservation laws such as the Euler and Navier-Stokes equations for problems including extreme shocks, shock-vortex interactions, and complex compressible turbulent flows.


Introduction
The use of discontinuous spectral element methods (DSEM) has grown in prevalence over the years due to their ability to achieve high-order accuracy while retaining geometric flexibility and a compact data structure suited for massively-parallel computing. As a result, DSEM offer many advantages for simulations of complex problems. However, their robustness is severely reduced for problems that exhibit discontinuities as the presence of spurious oscillations due to Gibbs phenomena can result in nonphysical solutions or the failure of the scheme altogether. Consequently, this lack of reliability is one of the limitations preventing the widespread adoption of these methods in the industry. To extend the use of DSEM to a wider variety of problems, various stabilization techniques have been proposed to increase the robustness of these schemes in the vicinity of discontinuities. A common goal of these shock capturing methods is to suppress numerical instabilities in the vicinity of a shock without degrading the accuracy of the underlying numerical scheme in regions where the solution is smooth.
These stabilization techniques can be broadly categorized as artificial viscosity, limiting, stencil modification, or filtering methods. The most ubiquitous approach is the addition of artificial viscosity to explicitly introduce numerical dissipation in the vicinity of a shock, a method pioneered by von Neumann and Richtmyer [1] with various advancements and alternative approaches over the decades [2][3][4][5][6]. Limiting methods have also shown promise [7,8], with the premise of these schemes generally relying on a combination of a constraint-satisfying low-order scheme and a constraint-violating high-order scheme. Additionally, modifying the numerical stencil to alleviate the issues of Gibbs phenomena has shown to be accurate and robust in other numerical settings [9][10][11][12], but the use of these methods is not very prevalent in the context of DSEM, primarily due to the incompatibility of the approach with the compact data structure, local compute, and geometric flexibility afforded by DSEM. In contrast, filtering [13][14][15], where spurious oscillations are removed by reducing high-frequency modes in the solution, is particularly attractive as a shock capturing method for DSEM as the filter can sharply resolve discontinuities with minimal computational cost and without sacrificing the computational efficiency of the underlying numerical scheme.
However, a typical drawback of many of these approaches is that they (1) do not necessarily guarantee that physical constraints on the solution are satisfied, which may lead to the failure of the scheme, (2) have free parameters which can require problem-and mesh-dependent tuning, a cost that cannot be afforded for large scale-resolving simulations, and/or (3) are not easily and efficiently implemented in the context of explicit DSEM on modern computing architectures. These issues motivate the development of shock capturing approaches for DSEM that can guarantee certain physical constraints are satisfied without requiring problem-dependent tunable parameters. In this work, we present an adaptive filtering approach for shock capturing in nodal DSEM to address these issues. By formulating physical constraints such as positivity and a local minimum entropy principle as constraints on the discrete solution, the filter strength is computed via a simple scalar optimization problem requiring only element-local information. Under some basic assumptions on the properties of the numerical scheme, the filtered solution is guaranteed to satisfy these constraints, resulting in an efficient and robust method for resolving discontinuous features without the use of problem-dependent tunable parameters. Furthermore, the proposed filtering approach does not appreciably degrade the efficiency and accuracy of the standard DSEM approach for smooth solutions.
The remainder of this paper is organized as follows. Preliminaries regarding the methodology are presented in Section 2. The proposed filtering approach is then presented in Section 3 followed by numerical implementation details in Section 4. The results of numerical experiments on the Euler and Navier-Stokes equations are then given in Section 5, followed by conclusions drawn in Section 6.

Preliminaries
To aid in presenting the proposed adaptive filtering approach, some preliminaries will be introduced in this section that briefly touch on the topics of hyperbolic conservation laws, minimum entropy principles, discontinuous spectral element methods, and modal filtering.

Hyperbolic Systems and Entropy Principles
The present work pertains to approximations of hyperbolic conservation laws (or hyperbolic components of mixed hyperbolic-parabolic conservation laws) of the form where ∈ ℝ is a vector-valued solution, is an arbitrary number of field variables, ( ) ∈ ℝ × is the flux, Ω ∈ ℝ is the domain, and is an arbitrary number of spatial dimensions. For simplicity, the domain is assumed to be periodic or the solution is compactly supported. If a solution of Eq. (1) is an entropy solution [16], an entropy inequality of the form can be posed, for which the inequality holds given any entropy-flux pair ( , ) [17] that satisfies = .
This notion of the entropy is generally referred to as a numerical entropy [18]. We utilize a formulation of Eq. (2) that is of opposite sign in comparison to the work of Tadmor [18] for consistency with a physical entropy. In smooth regions, the inequality is satisfied exactly (i.e., the entropy functional obeys a conservation law), whereas in the vicinity of shocks, the left-hand side attains a strictly positive value (i.e., an entropy source exists). For certain systems, entropy solutions of Eq. (1) satisfy a minimum principle on the entropy, i.e., for all Δ > 0 [18,19]. Assuming a finite propagation speed in hyperbolic systems, it is possible to restrict this condition over a local domain of influence to form a local minimum principle on the entropy. For a given point 0 , let 0 be some local domain of influence over the interval [ , + Δ ] -e.g., a closed -ball centered on 0 with a radius of max Δ , where max is some upper bound on the local maximum propagation speed of the system. A local minimum entropy principle can then be given as For numerical approximations of hyperbolic conservation laws, it can be advantageous to find schemes that can enforce this condition in scenarios where enforcing Eq. (2) is not feasible as it still guarantees some notion of physicality to the solution and can help alleviate instabilities in the numerical scheme [20].

Discontinuous Spectral Element Methods
For a nodal discontinuous spectral element approximation of Eq. (1) (e.g., discontinuous Galerkin [21], flux reconstruction [22], etc.), the domain Ω is partitioned into elements Ω such that Ω = ⋃ Ω and Ω ∩ Ω = ∅ for ≠ , as shown in Fig. 1. The approximate solution ( ) within each element Ω is given via a nodal approximation of the form where ∀ ∈ is a set of solution nodes, ( ) are their associated nodal basis functions that possess the property ( ) = , and is the set of nodal indices for the stencil. The shorthand notation is used to denote the discrete nodal value of the solution (i.e., = ( )). The order of the approximation, represented by ℙ for some order , is defined as the maximal order of ( ). The interaction between elements is communicated via the element interfaces Ω using a set of interface nodes ∈ Ω ∀ ∈ , where is a set of nodal indices for the interface stencil. In the work, we assume that the solution nodes are closed, such that the interface nodes are a subset of the solution nodes (i.e., ⊂ ) to avoid issues regarding interpolation for discontinuous solutions. Modifications to allow for open solution nodes are possible but are outside of the scope of this work. Given an interface point for some ∈ , let − denote the value of ( ) evaluated from the element of interest and let + denote the value of ( ) evaluated from the interface-adjacent element.
For an arbitrary element Ω of some discontinuous spectral element approximation, the divergence of the flux at some point can be approximated in terms of an interior flux component ( ⋅ Ω ( )) and an interface flux component For the interior component, the divergence of the flux is calculated through a collocation projection of the flux onto the solution nodes, given as for some discretization-dependent matrix . For the interface component, the divergence of the flux is calculated as a function of both the interior and exterior interface values, given as where is again some discretization-dependent matrix and ( − , + , ) is a common interface flux value dependent on the interior/exterior solution values ( − , + ) and their associated normal vector . This common interface flux is generally computed using exact or approximate Riemann solvers such as that of Rusanov [23] and Roe [24]. The semidiscretization of Eq. (1) can then be given as We assume that the spatial scheme is chosen such that where is the corresponding quadrature weight for the point and is the element-wise mean defined as for some arbitrary element Ω . This relation is recovered for nodal discontinuous Galerkin approximations with appropriate quadrature and flux reconstruction schemes utilizing the equivalent discontinuous Galerkin correction functions [25], the latter of which is the focus of this work. The ability of DSEM to preserve desirable properties on the elementwise mean is well documented in the literature due to its equivalency to first-order finite volume schemes [26][27][28].
Without presenting a comprehensive proof, we assume that properties such as positivity of certain convex functionals of the solution and a discrete local minimum entropy principle are satisfied by the element-wise mean if one utilizes explicit strong stability preserving time integration under some Courant-Friedrichs-Lewy (CFL) condition [29] with an appropriate choice for the Riemann solver. A more detailed description of these conditions and assumptions is presented in Section 4.

Modal Filtering
The approximation of the solution given by a nodal basis can be equivalently expressed by a modal expansion as where ( ) ∀ ∈ are a set of modal basis functions (e.g., Legendre polynomials, Koornwinder polynomials, etc.) and̂ are their corresponding modes. The modal basis is generally chosen such that the basis functions are orthogonal with respect to the inner product. With this formulation, a filtered solution can be defined as

Methodology
To enforce certain desirable properties of the systems in question (at least in a discrete sense), it is possible to formulate them as convex constraints on the solution. A common constraint for physical systems is that some convex functional Γ ( ) is non-negative across the domain (e.g., density and pressure for the Euler and Navier-Stokes equations, water height in the shallow water equations, etc.). This constraint is given by the condition where the choice of functional(s) is dependent on the system in question. However, simply enforcing positivity of these functionals is usually not enough to ensure a well-behaved solution in the vicinity of a discontinuity. This constraint must generally be accompanied by a more restrictive condition, such as some local discrete minimum entropy principle on the solution, given by the condition where ( ) is some convex entropy functional of the system in question and min is some local minimum entropy to be defined in Section 3.2. Enforcing this condition on the entropy tends to alleviate many of the issues regarding spurious oscillations in the vicinity of discontinuities [20], although this is not always guaranteed if this condition is not enforced for every possible entropy-flux pair.

Adaptive Filtering
To enforce these properties, we define an adaptive filtering operation with the goal of satisfying the constraints on the discrete filtered solution without the need for problem-dependent tunable parameters. The specific choice of filter is not particularly important as long as the filter meets the following criteria: • The filter is dependent on a single free parameter (i.e., ( ) = ( , )).
• There exists a minimum (or maximum) value of such that the filter recovers the unfiltered solution (i.e.,̃ ( ) = ( )).
• There exists a maximum (or minimum) value of such that the filter recovers the mean mode (i.e.,̃ ( ) = ).
The objective of the adaptive filter is to apply the minimum amount of filtering to the solution such that these constraints are met. We assume that for some arbitrary system, the discretization is chosen such that satisfies Eqs. (16) and (17), an assumption that is explored in more detail in Section 4 for the specific systems in this work. From this assumption, it can be seen that there exists at least one value of that recovers a filtered solution for which the constraints are satisfied exactly, i.e., Eqs. (16) and (17) are satisfied and at least one inequality becomes an equality. This value of can be calculated through simple root-finding methods such as the bisection algorithm, and, in practice, is generally unique, although this is not guaranteed. In this work, a second-order exponential filter [13] is chosen, given by the filter function where is defined as the maximal order of the modal basis function ( ). For this choice of filter, setting = 0 recovers the unfiltered solution and setting = ∞ recovers the mean mode, the latter of which can be approximated in a computational sense as =  (− log( )) for some value of machine precision . With this formulation, we define the filter strength using the minimum value of such that the solution abides by the positivity-preserving and discrete minimum entropy principle satisfying conditions, given by From a computational perspective, convergence to a local minima of is sufficient in the case that there exist multiple values of such that the constraints are satisfied exactly. A description of the approach for computing as well as the choice of functionals to enforce constraints upon for the various systems is given in Section 4. For sufficientlyresolved smooth solutions, the unfiltered solution is expected to already abide by these constraints [2], and therefore the standard DSEM approximation would be recovered. The proposed filtering approach is hereafter referred to as entropy filtering as the method effectively filters the modes of the solution that contribute to the violation of a minimum entropy principle.
Remark (Limiting). The proposed adaptive filtering operation can be considered somewhat similar to limiting-type approaches such as flux-corrected transport [7] and convex limiting [8]. In fact, for an unconventional filtering operation, given by the filter recovers the approach of Zhang and Shu [26], a linear convex limiting operation between a low-order solution ( = 0) and a high-order solution ( = 1). For a less trivial filter choice, this can be considered to be a nonlinear (and generally non-convex) limiting operation.

Entropy Constraints
Although positivity-preserving constraints are generally unequivocal for most physical systems, the notion of a minimum entropy principle is more ambiguous. Many systems are endowed with a multitude of entropy functionals ( ) [18] and the choice of min is not clearly defined. In this work, the minimum entropy principle is enforced on a numerical entropy that can be chosen arbitrarily from any entropy-flux pair ( , ) that satisfies Eq. (2) for the given system, the particular choice of which is posited to have a minor overall effect.
To calculate min , the discretization is assumed to be explicit in time under some standard CFL condition. With this assumption, the domain of influence of an arbitrary element Ω over a single temporal integration step can be considered to be strictly contained within the element and its direct Voronoi neighbors. Thus, min can be defined using information only from an element and its direct neighbors. Let * be defined as the minimum entropy within an element Ω , given as * = min where denotes the solution within the element Ω . Furthermore, let ( ) be the set of element indices which are face-adjacent with Ω , including Ω itself. The local minimum entropy min associated with the element Ω is then calculated as For elements adjacent to boundaries, the entropy of the boundary state is used as the adjacent entropy value. Assuming that these minimum entropy values are calculated prior to a temporal integration step and the filter is applied afterwards, this formulation of the entropy constraint enforces that the local discrete minimum entropy is non-decreasing in time across its domain of influence.

Governing Equations and Constraints
The efficacy of the entropy filtering approach was evaluated on hyperbolic and mixed hyperbolic-parabolic conservation laws. For the hyperbolic system, the compressible Euler equations were chosen, given in the form of Eq. (1) as where is the density, is the momentum, is the total energy, = ( −1) ( − 1 ∕2 ⋅ ) is the pressure, and = 1.4 is the ratio of specific heat capacities for air. The symbol denotes the identity matrix in ℝ × and = ∕ denotes the velocity. Positivity constraints were placed on the density and pressure, Γ 1 ( ) = and Γ 2 ( ) = , and the entropy functional was chosen as = log( − ), taken from the entropy-flux pair ( , ). Assuming an explicit strong stability preserving temporal integration scheme under some standard CFL condition with a solution that initially satisfies these constraints discretely, the element-wise mean at the next temporal step will satisfy these constraints if the interface fluxes are computed using an entropy-stable positivity-preserving Riemann solver [28,30] (e.g., Godunov methods [31], local Lax-Friedrichs flux [32], HLLC [33] with appropriate wavespeed estimates). These properties of the element-wise mean were shown for discontinuous Galerkin approximations in the works of Zhang and Shu [26], Zhang and Shu [27], Zhang et al. [28], and Chen and Shu [30].
Remark (Source terms). The extension of the proposed approach to hyperbolic systems with source terms is possible for discretizations that preserve the constraints on the element-wise mean. For positivity constraints, Zhang and Shu [34] showed these properties under a potentially more restrictive time step condition. For the entropy constraint, the contribution of the source term to the entropy over the temporal integration step would have to be evaluated to augment the min value, or if the entropy source is strictly positive, one may forego this modification and apply the proposed filter approach at the expense of a more-relaxed entropy constraint.
For the mixed hyperbolic-parabolic system, the compressible Navier-Stokes equations were chosen. As the assumption on the entropy of the element-wise mean is not necessarily satisfied with the mixed hyperbolic-parabolic discretization, the hyperbolic and parabolic components of the conservation law were isolated and treated separately using an explicit operator splitting approach [35]. For this system, the conservation law can be represented as where the solution and the inviscid (hyperbolic) and viscous (parabolic) components of the flux, denoted by the subscripts and , respectively, are given as where is the dynamic viscosity, = 0.73 is the Prandtl number for air, and = −1 ( − 1 ∕2 ⋅ ) is the specific internal energy. Identical positivity and entropy constraints are used for the Navier-Stokes equations as for the Euler equations. To enforce an entropy constraint for this system, at each temporal integration stage, the hyperbolic step is computed first, after which the entropy filter can be applied as with purely hyperbolic systems. The parabolic component of the temporal update is then added to the filtered solution. A final check is performed to ensure that the mixed hyperbolic-parabolic solution retains the positivity-preserving properties of the hyperbolic step. In the rare occasion that it does not, the filter is applied again using only positivity constraints. For both components of the flux, the boundary conditions were modified to ensure consistency with the system being solved (e.g., no-slip boundary conditions for the parabolic step were replaced with slip boundary conditions for the hyperbolic step).
Remark (Navier-Stokes equations). The notion of a minimum entropy principle is satisfied by the Navier-Stokes equations when considering the thermodynamic entropy (see Tadmor [18], Section 3). If the discretization for the parabolic component can be formed such that the minimum entropy principle is satisfied on the element-wise mean, the operator splitting approach can be neglected and the filter can be applied on the full hyperbolic-parabolic step with a significant reduction in computational cost. Alternatively, one may neglect the operator splitting approach without modifying the parabolic discretization to reduce the computational cost at the expense of the entropy constraint not necessarily being satisfied.

Discretization and Computational Framework
The proposed entropy filtering approach was implemented within PyFR [36], a high-order unstructured flux reconstruction [22] (FR) solver that can target multiple compute architectures including CPUs and GPUs. The FR framework can be cast in the form of Eq. (7) and Eq. (8) as and where ( ) corresponds to the Lagrange interpolating polynomial for the -th solution point and ( ) corresponds to the correction function for the -th interface point. The correction functions [37,38] posses the properties that where RT is the Raviart-Thomas space [39] of order . In this work, the correction functions are chosen such as to recover the nodal discontinuous Galerkin approach [21,22]. For a more in-depth overview and details on the extension to second-order PDEs, the reader is referred to Witherden et al. [40] and the references therein.
The solution nodes were distributed along the Gauss-Legendre-Lobatto quadrature points for tensor-product elements and the -optimized points [21] for simplex elements. Common interface flux values were computed using the HLLC Riemann solver [33] for the inviscid fluxes and the BR2 approach [41] for the viscous fluxes. Temporal integration was performed using a three-stage third-order strong stability preserving (SSP) Runge-Kutta scheme [42]. For the modal basis, orthogonal polynomials with respect to the unit measure were used (i.e., Legendre basis for tensor-product elements, Proriol-Koornwinder-Dubiner-Owens basis for triangles, etc.).

Filter Implementation
The implementation of the entropy filter was formulated as an element-wise scalar optimization problem. At each substage of the temporal integration method, a filtering operation was performed on the solution to enforce the positivity-preserving and minimum entropy constraints, the latter of which was computed using the solution at the previous substage. If the unfiltered solution satisfied the constraints, no filter was applied. Otherwise, the minimum necessary filter strength was calculated via 20 iterations of a bisection approach. Faster convergence could be obtained using more sophisticated root bracketing methods such as the Brent or Illinois methods, but these approaches were not explored in this work. For the density and pressure constraints, a minimum value of min = min = 10 −8 was enforced to ensure a non-vacuum state for the Riemann solver, such that the positivity constraints were instead implemented as Γ 1 ( ) = − min and Γ 2 ( ) = − min .
For the entropy constraint, a numerical tolerance of = 10 −4 was given, such that the constraint was instead implemented as This comparatively larger tolerance resulted in a slightly relaxed entropy constraint which was found to be beneficial for two reasons. First, the constraint was notably more prone to numerical precision issues due to the logarithm operation, particularly in the limit as → min or → min since small variations in these values could cause orders of magnitude more variation in the entropy. Secondly, marginally better resolution of flow features could be obtained by allowing slight undershoots in the entropy as strictly enforcing the entropy principle can degrade the accuracy of the solution [8,20]. A more detailed description of the computational implementation including pseudo-code is presented in Appendix A.

Results
The proposed entropy filtering approach was evaluated on a series of numerical experiments for the Euler and Navier-Stokes equations within a high-order flux reconstruction framework. For brevity, the solution for these systems is expressed in terms of a vector of primitive variables as = [ , , ] .

Sod Shock Tube
For an initial evaluation of the shock capturing capabilities of the proposed approach in the context of the Euler equations, the canonical case of the Sod shock tube was considered [43]. The problem assesses the ability of the approach in resolving the three main features of the Riemann problem: shock waves, rarefaction waves, and contact discontinuities. The domain is set to Ω = [0, 1] and the initial conditions are given as The density profiles at = 0.2 as predicted by a ℙ 3 and ℙ 5 FR approximation with 200 degrees of freedom are shown in Fig. 2. For both approximation orders, the results showed good agreement with the exact solution, with excellent resolution of the rarefaction wave and shock wave and minimal dissipation around the contact discontinuity. Furthermore, negligible spurious oscillations were observed in the vicinity of discontinuities. For a fixed number of degrees of freedom, marginally better results were obtained using the lower-order ℙ 3 approximation than the higherorder ℙ 5 approximation, particularly around the contact discontinuity. This effect can be attributed to the proportionally higher number of elements available for lower-order approximations, giving a more localized approach for the filter.
For a quantitative evaluation of the entropy filtering approach for discontinuous solutions, the convergence rates of the error against the exact solution were evaluated. For a given number of degrees of freedom , the point-mean 1 and 2 norm of the density error was defined as respectively. The convergence rates of the density error with respect to the number of elements for various approximation orders are shown in Table 1 and Table 2. The expected first-order convergence rate was generally obtained in both the 1 and 2 norm for all approximation orders, and the trend for the error was to decrease with increasing approximation order for a fixed number of elements. For a fixed number of degrees of freedom, shown in Table 3 and Table 4, marginally lower error was generally obtained with a lower-order approximation due to the previously mentioned effects.

Shu-Osher Problem
To assess the effects of the entropy filter for more complex problems including shock waves and smooth oscillatory behavior, the case of Shu and Osher [10]  The problem consists of a shock front propagating through a sinusoidally-perturbed density field, the interaction between which can induce instabilities in the flow field. However, these instabilities can be erroneously damped by overly dissipative shock capturing schemes. The predicted density profile at = 1.8 computed using a ℙ 3 FR approximation with 100 and 200 elements is shown in Fig. 3. A reference solution was computed using a highly-resolved exact Godunov-type solver [44]. The results show good resolution of the leading and trailing shock waves, and no spurious oscillations were observed. The instabilities in the field aft of the leading shock were also well-resolved, particularly with increasing resolution.

Isentropic Euler Vortex
For an extension to two-dimensional problems, the entropy filtering approach was initially applied to a smooth solution to verify the accuracy of the underlying DSEM was not detrimentally affected. The isentropic Euler vortex problem [11] was used as its analytic solution can be utilized to evaluate the convergence of the error. The initial conditions of the problem are given as with the parameters = 13.5 denoting the strength of the vortex, = 1.5 the radius, = 0, = 1 the advection velocities, and = 0.4 the freestream Mach number. The domain was set to Ω = [−10, 10] 2 , and a series of uniform quadrilateral meshes of size × was generated with periodic boundary conditions. After a single pass-through of the domain, the 2 norm of the density error was calculated as with the quadrature calculated on (2 ) 2 Gauss-Legendre nodes. This error is tabulated in Table 5 for a series of experiments with varying mesh resolution and approximation order. The convergence rate of the error was generally in the range of to + 1 for the varying approximation orders, on par with the theoretical rate of + 1. These findings indicate that the entropy filter does not appreciably degrade the accuracy of the underlying DSEM for smooth solutions.  To verify the computational efficiency of the entropy filter for smooth solutions where the unfiltered solution remains stable, the compute time for a solution without filtering and with a varying number of filter iterations was compared. The comparison was performed on an NVIDIA V100 GPU over 5 flow-throughs of the domain on an = 40 mesh with a ℙ 3 approximation and Δ = 2⋅10 −4 . For 5, 10, and 20 iterations of the filter, the relative computational cost increase (as measured by relative time-to-solution) was 1.1%, 1.7%, and 2.0%, respectively. These findings indicate that the entropy filter does not have a notable detrimental impact on the computational efficiency of the scheme for smooth solutions as it is primarily inactive. Furthermore, it also indicates that further improvements to the efficiency are possible through more sophisticated iterative solvers that require fewer iterations.

Double Mach Reflection
The double mach reflection problem of Woodward and Colella [45] was subsequently used to evaluate the ability of the entropy filter to resolve strong discontinuities in multiple dimensions. This case consists of a Mach 10 shock impinging on a 30 degree ramp and results in multiple strong shock-shock and shock-contact interactions. The problem is solved on the domain Ω = [0, 4] × [0, 1] with the initial conditions At the left boundary and the bottom boundary for < 1∕6, the solution was set to the post-shock state . No-slip adiabatic wall boundary conditions were applied for the bottom boundary for ⩾ 1 6 . At the right boundary, the solution was set to the pre-shock state . For the top boundary, the exact solution is enforced, given as The contours of density at = 0.2 as predicted by a ℙ 3 FR scheme on a 2400 × 600 mesh are shown in Fig. 4. The results show sub-element resolution of discontinuities without the presence of spurious oscillations. Furthermore, the application of the filter did not excessively dissipate the Kelvin-Helmholtz instabilities along the contact line, indicating that the filter is not erroneously dissipating small-scale features. To verify this, the distribution of the filter parameter is also shown in Fig. 4 overlaid on the isocontours of density. The filter was primarily active in the leading shock fronts, with minimal activation within the small-scale structures along contact line. Furthermore, even along the shock front, a relatively small value of was required, significantly less than the value corresponding to the recovery of the mean mode.

Kelvin-Helmholtz Instability
To further evaluate the ability of the entropy filtering approach to resolve small-scale flow features in the vicinity of discontinuities, the roll-up of a Kelvin-Helmholtz instability was simulated. The problem is solved on the domain Ω = [−0.5, 0.5] 2 with the initial conditions No explicit seeding of the instabilities was performed -the instabilities originate from roundoff errors in the solver. The contours of density at = 2 as predicted by a ℙ 4 FR scheme with various mesh resolutions are shown in Fig. 5.
The roll-up of the vortices was well-resolved by the approach, with complex small-scale vortical structures beginning to appear with increasing mesh resolution. No spurious oscillations were observed in the vicinity of the discontinuities, and more subtle features such as pressure waves were not excessively dissipated.

Mach 800 Astrophysical Jet
As a verification of the positivity-preserving properties of the filter for very extreme conditions, the case of highspeed astrophysical jets is considered. The test case, introduced by Balsara [46], consists of a Mach 800 jet in the presence of an ambient gas. The problem setup is identical to the work of Wu and Shu [47] in which a provablypositive third-order discontinuous Galerkin approach is used, although the magnetic field is neglected in this case. A half domain Ω = [0, 0.5] × [0, 1.5] is considered with symmetry (slip adiabatic) boundary conditions along the -axis. The domain is filled with an ambient gas of density 0.1 , zero velocity, and unit pressure. For the = 0 boundary, the inlet region is defined on ⩽ 0.05, and the solution is set to = [ , 0, 800, 1] which yields a Mach number of 800 with respect to the inflow gas. The remaining boundary conditions are set as free.
The contours of density at = 0.002 as predicted by a ℙ 3 FR scheme with a coarse (200×600) and fine (800×2400) mesh are shown in Fig. 6. Excellent resolution of the leading shock wave was obtained, and the small-scale structures in the vicinity of the cocoon/jet interface were well-resolved even with the coarse mesh. Additionally, the distribution of the filter parameter is shown in Fig. 6 overlaid on the isocontours of density. For both the coarse and fine mesh, the filter was primarily active at the leading shock region. Some regions of activation within the cocoon/jet interface were observed for the coarse mesh, but this behavior was reduced with increasing resolution, such that minimal activation away from the leading shock was observed with the fine mesh.
Given a similar case setup, a comparison can be made between the proposed approach on the coarse mesh and the mildly-magnetized results of Wu and Shu [47] (Fig. 6.i), which were obtained by the linear limiting approach of Zhang and Shu [26]. Significantly better resolution was obtained using the entropy filter, particularly with regards to the small-scale features near the cocoon/jet interface, although some of this may be attributed to the stabilizing effect of the magnetic field in their case. However, this suggests that the nonlinear limiting performed by the entropy filter offers noticeable advantages in comparison to a linear limiting approach with minimal computational overhead. Contours are reflected about the -axis. View is reoriented such that the + direction is shown left-to-right.

Taylor-Green Vortex
The entropy filtering approach was then extended to turbulent compressible flows through the Navier-Stokes equations. To verify that the filter does not unnecessarily dissipate small-scale turbulent fluctuations in the absence of shocks, the proposed approach was applied to the subsonic Taylor-Green vortex at a Reynolds number of 1600, a canonical fluid dynamics problem for studying vortex dynamics and turbulent transition and decay [48]. The problem is solved on the periodic domain Ω = [− , ] 3 with the initial conditions  (2 )) (cos(2 + 2)) where 0 = 1∕ 2 for a reference Mach number = 0.08. Given the unit density and velocity, the dynamic viscosity is set to 1∕1600 to recover a Reynolds number of 1600.
The quantity of interest for this problem is the dissipation rate of the kinetic energy in the flow. The non-dimensional integrated kinetic energy can be defined as where = 8 3 is the volume of the domain. From this, a dissipation rate based on the kinetic energy can be calculated as A similar measure of the dissipation can be obtained through a scaled form of the non-dimensional integrated enstrophy, defined as where is the vorticity and = 2 is the scaling factor. For purely incompressible flows, these two quantities are equal, but for compressible flows, they differ by the contribution of the deviatoric strain and pressure dilatation to the dissipation. At low Mach numbers, the enstrophy-based dissipation can reasonably approximate the kinetic energybased dissipation for well-resolved flows. The prediction of these two quantities as computed by the entropy filtering approach with varying resolution and approximation order is shown in Fig. 7 in comparison to the DNS results of van Rees et al. [49]. For visualization of the kinetic energy-based dissipation, a moving-average smoothing operation was performed prior to computing the temporal derivative to reduce oscillations. For a ℙ 3 FR approximation, the number of degrees of freedom was varied from 96 3 -160 3 . Relatively good agreement was observed between the kinetic energy-based dissipation and the reference data across this entire range of resolution, with the most resolved case showing negligible deviation from the reference. For the enstrophy-based dissipation, the profiles evidently showed convergence to the reference with increasing resolution, and for a given resolution, showed similar results to unfiltered approaches [50]. When fixing the degrees of freedom to ∼120 and varying the approximation order, improvements in the prediction of both the kinetic energy-and enstrophy-based dissipation were observed with increasing approximation order. Additionally, nearly identical results were obtained without filtering (see Trojak et al. [51], Fig. 3b), which supports the presumption that the filter is predominantly inactive for turbulent flows without discontinuities.

Viscous Shock Tube
To assess the efficacy of the entropy filtering approach for predicting shock-boundary layer interactions, a twodimensional viscous shock tube was simulated. The problem, introduced by Daru and Tenaud [52], is solved on the half-domain Ω = [0, 1] × [0, 0.5] with the initial conditions For the top wall, a slip adiabatic wall boundary condition is applied to enforce symmetry. No slip adiabatic wall conditions are applied for the remaining walls. This case contains the standard features of the Riemann problem, namely a rarefaction wave, a contact discontinuity, and a shock wave, in addition to viscous interactions between these features and the wall. This interaction forms a complex lambda shock impinging on a viscous boundary layer, and as a result, presents a suitable test case for shock-boundary layer interactions. Numerical Schlieren results at = 1 computed using a ℙ 4 FR approximation on a uniform 800 × 400 mesh are shown in Fig. 8 for varying values of the dynamic viscosity . As the Reynolds number was increased, progressively smaller-scale features became evident in the flow, particularly along the contact line and in the separation region. The discontinuities as well as the small vortical structures in the flow were both well-resolved. The predicted results show relatively good agreement with the results of Guermond et al. [53], although some small-scale features at higher Reynolds numbers were notably more pronounced in the present work. These results suggest that the entropy filter used in conjunction with an operator splitting approach is an effective tool for approximating mixed hyperbolic-parabolic systems that exhibit discontinuities as well as complex small-scale flow structures.
To investigate the behavior of the filter for viscous flows, the distribution of the filter parameter overlaid on isocontours of density is shown in Fig. 9 for the case presented in Fig. 8d. The distribution for the present case was more sporadic than in the inviscid cases, with more frequent activation near shock-boundary layer interactions (some of which are more easily observed with the Schlieren diagram in Fig. 8d than with the isocontours in Fig. 9). It was also observed that for a small portion of the elements within the mesh, the addition of the viscous component caused the initially positivity-preserving solution to violate the positivity constraints. Approximately 16 elements required additional filtering, but the amount of filtering required was very minor, with values of on the order of 10 −4 .

Transonic Delta Wing
As a final evaluation of the proposed approach for complex aeronautical applications including three-dimensional high Reynolds number flows computed on unstructured meshes, the test case of large eddy simulation around a transonic VFE-2 delta wing was considered. The geometry, introduced by Chu and Luckring [54], consists of a sharp leading edge delta wing with a sweep angle of 65 • and a thickness to root chord ratio of 3.4%. To match the experimental setup of Konrath et al. [55], the Reynolds number is set to 3⋅10 6 based on a unit root chord (i.e., 2⋅10 6 based on the mean aerodynamic chord), and the freestream Mach number is set to 0.8, yielding locally supersonic turbulent flow on the suction side of the wing. For these geometry and flow conditions, the separation and roll-up of the primary vortex is triggered by the sharp leading edge which immediately transitions into turbulence. Secondary and tertiary vortices are then subsequently formed through near-wall viscous interactions. An angle of attack of 20.5 • was chosen as these conditions yield strongly nonlinear behavior on the suction side as the primary vortex is on the precipice of vortex burst.
A second-order unstructured mesh was generated by extruding a triangular surface mesh in the near wall region and filling the remaining domain with tetrahedral elements. The resolution of the surface mesh was approximately uniform with an average edge length of ℎ = 0.003, and the near-wall prismatic region was extruded over a length of 0.02. To simplify the meshing, the sting and fairing from the experimental model were neglected. Since the separation and roll-up of the primary vortex is primarily governed by inviscid effects due to the sharp leading edge [56], the wall-normal resolution requirements are relaxed in comparison to standard wall-resolved large eddy simulation [57]. As such, the wall-normal resolution was set such that the + value of the first solution point away from the wall was approximately 30 (i.e., the element + was ∼120) based on equivalent flat plate conditions matching the freestream flow and root chord. In the separated region above the suction side, the resolution of the tetrahedral mesh was such that the average edge length was approximately ℎ = 0.01. The overall mesh consisted of 4.1 million elements, yielding a total of 104 million degrees of freedom per component with a ℙ 3 approximation. At the wing, no slip adiabatic wall boundary conditions were applied, and characteristic Riemann invariant boundary conditions were used for the farfield.
At startup, the approximation order and Reynolds number were progressively increased until the initial transients were convected away. Averaging was then performed for a time period corresponding to five flows over root chord. The results of a ℙ 3 FR approximation are presented in Fig. 10, showing instantaneous Q-criterion isosurfaces colored by local Mach number and time-averaged surface pressure coefficient contours. The Q-criterion isosurfaces showed the canonical vortex structure of delta wing flows with the primary and secondary vortices evident in the flow. Due to the relatively high angle of attack, a local Mach number of approximately 2 was observed near the primary vortex core, and indication of the onset of vortex burst was seen near the trailing edge. Furthermore, the average surface pressure coefficient contours clearly showed the presence of a primary and secondary vortex, and good agreement was observed between the predicted results and the experimental results of Konrath et al. [55].
The surface pressure coefficient at various streamwise locations is presented in Fig. 11 in comparison to the experimental data. At ∕ = 0.2, 0.4, and 0.6, excellent agreement was observed, both in terms of the location and the magnitude of the primary and secondary vortex pressure peaks. At ∕ = 0.8, the location and magnitude of the primary vortex pressure peak was well-resolved, but a slight overprediction in the magnitude of the secondary vortex Figure 10: Isosurface of instantaneous Q-criterion colored by local Mach number (left) and time-averaged surface pressure coefficient contours (right) for the transonic VFE-2 delta wing computed using a ℙ 3 FR approximation.
pressure peak was seen. This effect may be attributed to the strong sensitivity of the flow regime prior to vortex burst in conjunction with the slight discrepancies between the experimental and computational geometry at the trailing edge of the wing. Overall, the computational results were in-line with the experimental data, indicating that the entropy filtering approach can be successfully used for complex supersonic turbulent flows.  Figure 11: Surface pressure coefficient at various streamwise locations for the transonic VFE-2 delta wing computed using a ℙ 3 FR approximation. Experimental results of Konrath et al. [55] shown for reference. Dotted line denotes the sonic pressure coefficient. Spanwise extent is normalized by the local semispan .

Conclusions
In this work, a novel adaptive filtering approach for shock capturing in discontinuous spectral element methods is presented. The proposed entropy filtering approach formulates the filtering operation as an element-wise scalar optimization problem designed to enforce constraints such as positivity and a local discrete minimum entropy principle. As a result, this method is free of problem-dependent tunable parameters and can be easily implemented on arbitrary unstructured meshes with minimal computational effort. To evaluate the efficacy of the proposed filtering approach, a series of numerical experiments was performed on multi-dimensional hyperbolic and mixed hyperbolic-parabolic conservation laws such as the Euler and Navier-Stokes equations. Even for extremely strong shocks, sub-element resolution of the discontinuities was generally obtained with minimal spurious oscillations, and high-order accuracy was recovered for smooth solutions. Furthermore, for more complex problems including shock-vortex interactions and compressible turbulent flows, the filter was able to robustly resolve discontinuities without unnecessarily dissipating small-scale vortical structures in the flow. These results indicate that the proposed entropy filtering approach can be an effective and robust technique for shock capturing in discontinuous spectral element approximations of hyperbolic and mixed hyperbolic-parabolic conservation laws. In future work, the proposed approach can be improved by optimizing the calculation of the filter strength through more sophisticated root bracketing methods and stopping criteria, increasing the locality of the filtering kernel to mimic that of methods such as artificial viscosity and subcell approaches, modifying the approach to remove the need for operator splitting in mixed systems, and applying a more robust method for calculating the entropy tolerance.