Numerical Investigation of the Inviscid Taylor-Green Vortex Using an Adaptive Filtering Method for a Modal Discontinuous Galerkin Method

Implicit Large Eddy Simulation and under-resolved Direct Numerical Simulation bypass the complexity and uncertainty of turbulence modelling by using the numerical dissipation of the scheme as a subgrid scale model. High-order methods allow for more accurate capturing of smaller scale structures but suffer from energy pile-up in the higher modes which leads to instability in under-resolved applications. This work presents a filtered modal Discontinuous Galerkin method which adaptively determines the filter strength, avoiding unnecessary degradation of accuracy while maintaining stability. The method is applied to the inviscid Taylor-Green Vortex, a challenging test case which exhibits under-resolved turbulence for which few published results exist. This goal of this work is to present the adaptive filtering method which achieves robustness and accuracy despite a low number of degrees of freedom, as well as to publish a quantity of relevant data for the inviscid TGV problem.


Introduction
The wide range of spatial and temporal scales present in turbulence makes it unfeasible to fully simulate, also known as direct numerical simulation (DNS).Current approaches to computing turbulent flows include averaging the governing equations (RANS), which requires some type of modelling to provide closure (k-omega, k-epsilon, etc.); and Large Eddy Simulation (LES), which uses a subgrid-scale (SGS) model to model the smaller scales while solving a filtered set of equations.RANS is unsuitable for unsteady calculations, and the design of models for both LES and RANS presents nontrivial difficulties.The most popular SGS model currently, due to its simplicity, is the Smagorinsky-Lilly model (Smagorinsky 1963), which calculates the turbulent viscosity based only on a tunable constant, the filter width, and the large-scale strain rate.The constant must be adjusted based on the flow conditions and numerical scheme, and the model's largest shortcoming is that it is dissipative even for well-resolved flows.An improvement to the model is the dynamic Smagorinsky model (Germano et al. 1991), which adjusts the tunable constant based on the subgrid and CONTACT Dean Yuan dean0927@gmail.comsubfilter stresses.However, it can suffer from stability issues (Ghosal et al. 1995).Other models include the Spectral Vanishing Viscosity (Karamanos and Karniadakis 2000;Tadmor 1989), which adds a viscosity term whose amplitude decreases as the wavenumber decreases and has seen relative success, and scale similarity models (Bardina, Ferziger, and Reynolds 1980) which assume similarity between the subgrid-scale and resolved scales, but require additional dissipation models to maintain stability.These models have also been shown to produce very good results, but are more computationally costly.Other approaches include the variational multiscale (VMS) (Hughes, Oberai, and Mazzei 2001) which restricts the influence of the SGS to the small scales, and the Wall-Adapting Local Eddyviscosity model (WALE) (Nicoud and Ducros 1999) which achieves accurate scaling at the wall.
For an LES method to be successful, a numerical method which has very low amounts of dissipation and dispersion is required.The discontinuous Galerkin (DG) method has become a popular choice due to its high accuracy and parallelizability, and has been used for LES in a handful of papers.In Van Der Bos and Geurts (2010), a parametric study is done on the Smagorinsky parameter and it was found that the best results are achieved when the DG order is as high as possible and numerical dissipation from the Riemann solver is as low as possible.In Beck et al. (2016), it was found that aliasing errors dominates the tuning of the constants for Smagorinsky and VMS, and that implicit models showed better agreement with DNS results.In Chapelier, de la Llave Plata, and Lamballais (2016), a VMS method was developed for a modal DG method with over-integration and results are presented for the inviscid TGV at 48 3 resolution as well as the viscous TGV case.In de la Llave Plata, Couaillier, and Pape (2018), a modal DG method was used to study wall-bounded turbulent flows using uDNS and WALE, and found that WALE was only slightly more accurate than the uDNS for the 2D periodic hill problem, but an over-integration approach alone was not enough to stabilise the uDNS method at Re 19,000.
Recently gaining popularity is the implicit LES (ILES) approach, also known as embedded or implicit turbulence modelling (Drikakis 2002) or Monotone Integrated Large Eddy Simulation (MILES) (Boris et al. 1992), in which the inherent dissipation of the numerical method plays the role of the SGS model.This approach has the benefit of lower complexity due to the lack of an SGS model while also removing the uncertainty involved with modelling and its interaction with numerical errors, and has been applied to the DG method in a variety of publications (Beck et al. 2014;Carton de Wiart and Hillewaert 2015;Carton de Wiart et al. 2015;Fernandez, Nguyen, and Peraire 2017;Frère et al. 2015;Murman et al. 2016;Renac et al. 2015;Uranga et al. 2009Uranga et al. , 2011)).In Moura et al. (2016), a comparison is made of the iLES-DG methods to the LES method based on projection onto local basis functions (LES-PLB) Pope (2001), where the resolved scales are represented by a Galerkin projection and the residual represents the unresolved scales, which require modelling.Moura et al. (2016) shows that the DG formulation can be written as an LES-PLB method with terms representing implicit subgrid-scale modelling.However, no conclusions are able to be drawn about the suitability of DG as an iLES method based on the formulation alone.
Regardless of whether explicit or implicit modelling is employed, the approximation error of the numerical scheme plays an important role in accurately simulating the effects of turbulence.The use of high order methods such as continuous or discontinuous Galerkin (CG/DG), spectral volume, spectral difference, or flux reconstruction allows for more accurate capturing of the small scale turbulent structures due to their low dissipation/dispersion properties.This paper focuses on the discontinuous Galerkin method with implicit modelling, also known as underresolved DNS, which has gained a lot of popularity in recent years.However, the DG method suffers from aliasing instabilities, a well known issue caused by integration errors, as existing quadrature rules are unable to exactly integrate the nonlinear fluxes of the Euler/Navier-Stokes equations.Several approaches to address this issue are present in the literature, such as consistent or overintegration, i.e. increasing the number of quadrature points to decrease the integration error (Moura et al. 2017); filtering, which uses a filtering mechanism to dissipate out instabilities (Flad, Beck, and Munz 2016;Gassner and Beck 2013;Hamedi and Vermeire 2022); and the use of split form or entropy stable formulations (Chan et al. 2022;Gassner, Winters, and Kopriva 2016;Winters et al. 2018).Over-integration increases drastically the computational cost of the method; as many as 2-4 times the number of points are required (Moura et al. 2017), while still not guaranteeing stability.Filtering indiscriminately leads to excessive dissipation, increasing the error of the scheme.Finally, entropy stable formulations are greatly promising although they present a nontrivial amount of implementation difficulty.In this work we present an adaptive filtering approach, and follows a methodology similar to that presented as the Adaptive Dissipation/Dispersion Adjustment (ADDA) method (Fernandez-Fidalgo et al. 2021;Nogueira et al. 2020), originally developed for FV methods.Other filtering based approaches can be found in Gassner and Beck (2013), which presents results for a DGSEM applied to the Taylor-Green Vortex; Flad, Beck, and Munz (2016), which uses a energy transfer balance criterion to determine when to apply a modal cutoff filter; Hamedi and Vermeire (2022), which evaluate various exponential filter strengths and finds that those resembling cutoff filters perform the best; and Dzanic, Trojak, and Witherden (2023), which uses an entropy criterion to iteratively determine the filter strength.
Nodal schemes are most common for DG implementations, where an interpolating polynomial such as the Lagrange is used to represent the approximate solution within each element and the degrees of freedom are the values at various nodes within each element.This approach has the benefit of computational efficiency when the nodes are chosen to collocate with the Gaussian integration points; however, the method is no longer collocated when over-integration is utilised.We choose instead a modal approach, where the degrees of freedom are the coefficients of a Taylor basis.This approach has the benefit of not being tied to the choice of basis or quadrature rule, which allows for arbitrary element shapes and nonconforming polynomial orders with no additional implementation.We find that we achieve relatively low amounts of dissipation with a much smaller number of degrees of freedom, and achieve robustness with an Adaptively Filtered Discontinuous Galerkin (AF-DG) method.
All the schemes are developed in the open source UCNS3D solver (Antoniadis et al. 2022;UCNS3D 2022), and we assess their performance using the inviscid Taylor-Green vortex problem.The paper is organised as follows: in Section 3 we introduce the numerical framework used to describe the AF-DG framework.The numerical results are presented in Section 4. Section 5 describes the conclusions drawn from this study.

Governing Equations
We examine the 3-dimensional Euler hyperbolic conservation laws for inviscid, compressible flow: where U = U(x, t) is the vector of conserved variables, x = (x, y, z) denotes the coordinates of a point of the domain , and is the non-linear flux tensor: The system of equations is closed with the ideal gas equation of state:

Discontinuous Galerkin Method
In DG methods, higher orders of accuracy are attained by means of a high-order polynomial representation of the local element solution.Similar to FV schemes, as the solution is non-unique at the element interfaces, the flux terms are computed using an approximate Riemann solver.
Consider the unsteady non-linear hyperbolic system of conservation laws on a 3D domain , written in its conservative form, as in Equation (1).
The physical domain consists of any combination of conforming tetrahedra, hexahedra, prisms, or pyramids in 3D, and quadrilaterals or triangles in 2D.All the elements are indexed by a unique mono-index i.The global solution is discretely approximated for each cell i by a set of locally continuous piecewise solutions U i (x).The local approximation of the solution within each cell i, is given by a linear combination of a basis functions φ k i (x) as follows: The vector U k h,i (t) denotes the unknown degrees of freedom of the solution that are advanced in time.We take the weak form of the DG formulation, which is obtained multiplying Equation (1) by a smooth test function ψ(x) integrating over each element V i ∈ and performing an integration by parts: where S i is the boundary of V i and the non-uniquely defined flux F(U i (x)) on each surface S i is replaced by a numerical flux function , n are the solution values from the internal and external sides of the surface and the normal vector, respectively.Adopting a smooth test function from the same space of the basis function, thus , the DG formulation for each element i is given by: A popular choice of scheme is the nodal collocation scheme, where the DoFs are the solution values at various nodes within each element.In this work a modal formulation is used, which allows for the same basis to be used for any element shape.In the modal formulation, the unknowns to be solved are the polynomial expansion coefficients, and we choose for the polynomial basis functions the cell-centred Taylor series expansion (Boscheri and Dimarco 2022;Gaburro et al. 2020;Luo, Baum, and Löhner 2008;Luo et al. 2010;Maltsev, Skote, and Tsoutsanis 2024): where 0 ≤ p k + q k + r k ≤ N. The index i ranges over the total number of elements, and x c i , y c i , z c i are the coordinates of the cell centre of the considered element i.The basis functions are intentionally designed such that the higher-order modes vanish when integrating within the considered cell i, which is convenient when pursuing conservation and seamless interaction between DG and FV.Of great importance is the normalisation of the basis function with a coefficient h, in this work taken as the square root of target cell's volume, used to improve the mass matrix conditioning.This basis is hierarchical but non-orthogonal.
The number of DOFs K depends upon the degree of the polynomial N and the spatial dimension d.For a tensor product nodal collocation scheme using the Gauss-Legendre (GL) or Legendre Gauss-Lobatto (LGL) nodes on the reference hexahedron, the number of DOFs, and equally, the number of integration points, is LGL schemes have an error in the mass matrix (Kopriva and Gassner 2010).For modal bases on arbitrary shapes, the number of DOFs can be calculated from the following expression: Table 1 lists the number of DoFs and quadrature points for a hexahedra for modal and nodal methods, with the nodal quadrature points defined by the tensor product GL points and the modal ones defined by Williams, Shunn, and Jameson (2014), showing that a much lower number of DoFs and similar number of quadrature points are used in modal formulations.For example, for a 5th-order polynomial, the number of DoFs is 56 for a modal basis and 216 for nodal.
Another benefit of the modal formulation is the ability to handle arbitrary unstructured meshes and element shapes without utilising different basis functions or quadrature rules for each element type.The integration of the basis functions in physical space uses a decomposition of elements into triangles or tetrahedra.This strategy has been previously documented in other studies (Boscheri and Dimarco 2022;Gaburro et al. 2020;Titarev, Tsoutsanis, and Drikakis 2010;Tsoutsanis, Titarev, and Drikakis 2011) and differs from the approaches utilised in other methods, such as tensor product Legendre-Gauss or Legendre-Gauss-Lobatto rules which are highly popular for nodal schemes.In our implementation, hexahedra are decomposed into 6 tetrahedra, and a quadrature rule of order N (Williams, Shunn, and Jameson 2014) can be used on each decomposed element, while rules of order N + 1 must be used for triangles, quadrilaterals, and tetrahedra.Additional element shapes have not yet been tested, but are expected to behave similarly.The mass matrix for each cell is given by: The non-uniquely defined intercell fluxes in Equation ( 6) are determined in the same manner as in the FV framework, i.e. with an exact or approximate Riemann solver (Cockburn and Shu 1989).The semi-discrete system for each element i in the domain can be written as: where Fn im is the numerical flux normal to the interface between cell i and each neighbouring cell m, N f is the number of faces comprising the boundary of each cell i, N qp is the number of face quadrature points, and |S s | is the surface area of the each decomposed face s ∈ N s .N v is the number of the volume Gaussian quadrature points for each decomposed element d ∈ N d , and x i,g and ω g are the numerical integration points and weights as specified in Williams, Shunn, and Jameson (2014), with the points being given as the coefficients of a linear combination of the vertices of each decomposed simplex, and |V d | is the volume of each decomposed simplex.The surface integral is treated similarly in 3 dimensions.U n im,L (x im,α , t) and U n im,R (x im,α , t) are the discontinuous solutions for cell i and each neighbouring cell m at the prescribed face Gaussian quadrature points respectively; while x im,α and ω α correspond to different face Gaussian integration points and weights respectively.Equation ( 10) can be written as: where U i is the approximate solution, M i is the mass matrix and R i (U i ) is the residual vector.

Fluxes & Time Advancement
The flux approximation and temporal discretization techniques are shared between the FV and DG methods.For the inviscid fluxes, the approximate Riemann solver of Rusanov (1962), also known as the local Lax-Friedrichs, is employed unless otherwise stated.Moura, Sherwin, and Peiro (2015) reports better performance and robustness from the HLLC and Roe solvers; however, as only the Rusanov solver is used in many of the publications with TGV results (Chapelier, De La Llave Plata, and Renac 2012;Gassner and Beck 2013;Hamedi and Vermeire 2022), this work does as well for ease of comparison.The solution is advanced in time by a third-order TVD Runge-Kutta method.Because the mass matrix for the Taylor basis employed is non-diagonal, special care must be taken when filtering, as the degrees of freedom are coupled.We employ the procedure detailed in Kuzmin (2013), applying the filtering to the time derivative as well as the solution: where Ũ(n+k) i is defined as where is the filtering operation, M C is the consistent mass matrix as defined in 9, M L is the lumped mass matrix i.e. the diagonal elements of M C , and the time step t is selected according to: where S i is an estimate of the maximum in absolute value of the propagation speed in the cell i, h i is a characteristic length of the element i, and CFL refers to the Courant-Friedrichs-Lewy condition.
All the schemes developed are implemented in the UCNS3D CFD code (Antoniadis et al. 2022) which is written in object-oriented Fortran 2003, employing MPI message passing interface (MPI), and the Open Multi-Processing (OpenMP) application programming interface (API), and the reader is referred to Tsoutsanis, Antoniadis, and Jenkins (2018) and Tsoutsanis (2017) for more details on implementation and performance benchmarks.

Adaptive Filtering
As previously mentioned, it is well documented that methods which use a weak form of the equations such as DG suffer from so-called 'aliasing instabilities' due to the inexact integration of nonlinear terms.In the Euler/Navier-Stokes equations, this stems from the convective flux.In under-resolved settings, this leads to an unphysical 'build-up' of energy in the high modes, where it is trapped and unable to dissipate.The standard approach to alleviating this issue is either through 'over-integration', i.e. increasing the quadrature precision, or through filtering of the high modes, thus injecting artificial dissipation.Over-integration is computationally expensive and does not guarantee stability (Moura et al. 2017) and naive filtering is overly dissipative; therefore, various adaptive filtering approaches have been proposed (Dzanic, Trojak, and Witherden 2023;Fernandez-Fidalgo et al. 2021;Flad, Beck, and Munz 2016;Nogueira et al. 2020).
In Moura et al. (2017), the inviscid TGV case is run with an over-integrated nodal DG method, and it was found that for higher-order discretizations (N>5 or 6, depending on the Riemann solver), the simulation was unstable even when the number of quadrature points was increased to 4 times the number of DoFs.The reason for this instability is not fully understood, and it is hypothesised in Moura et al. (2017) that the sharper dissipation characteristics of higher-order DG discretizations creates a larger build-up of energy at the dissipation range due to the bottleneck effect (Coantic and Lasserre 1994;Falkovich 1994).
In Dzanic and Witherden (2022) and Dzanic, Trojak, and Witherden (2023), an exponential filter where the filter strength is determined iteratively in order to satisfy positive pressure and density as well as a minimum entropy constraint.The approach presented in this paper uses a similar procedure, but uses the energy ratio rather than entropy as the constraint In addition, a linear filter is used as no benefit to using an exponential filter has been shown in the literature: Hamedi and Vermeire (2022) performs a study on optimal exponential filter parameters by iteratively running an inviscid TGV case with a non-adaptive, globally filtered DG method and finds that the 'best' filters, defined as the weakest filter which maintains stability up to t = 20, resemble cutoff filters which filter the highest 1/3 or 1/2 modes, depending on the polynomial order.However, no kinetic energy, dissipation rate, enstrophy, or other results for the inviscid TGV are shown.
This work uses the ratio of high mode to low mode energy as the criterion for determining when to filter the solution, as in Li and Tsubokura (2017), Tantikul and Domaradzki (2010), Nogueira et al. (2020) and Fernandez-Fidalgo et al. (2021).In Tantikul and Domaradzki (2010), the incompressible Navier-Stokes equations are solved for turbulent channel flow, and filtering is activated when the energy ratio exceeds some theoretical amount.In Nogueira et al. (2020), the energy ratio is calculated using two different widths of a Moving Least Squares based filter and adjustment is made via the dissipative term of the Riemann solvers of Roe or Rusanov in a finite volume method for the compressible Navier-Stokes equations.In Fernandez-Fidalgo et al. (2021), the dissipative term of a finite difference WENO scheme for the incompressible Navier-Stokes equations is adjusted using the same energy ratio calculation.In our experiments, adjusting the dissipation of the Riemann solver was not enough to stabilise the simulations, which are Euler rather than Navier-Stokes.In Flad, Beck, and Munz (2016), the rate of change of the high mode and low mode energy is considered for the viscous TGV case; however, in our experiments, this energy transfer balance indicator was not effective for stabilising the inviscid TGV case.
In our approach, we iteratively determine the filter strength based on the ratio of high mode to low mode energy, where the solution is not filtered at all if the ratio is beneath a certain threshold and filtered down to a first-order polynomial at maximum dissipation.In practice, very few cells are marked for filtering and the filter strength is generally weak.
In order to determine the filter strength, we calculate within each element at each Runge-Kutta stage the ratio of high mode to low mode energy, where the high mode energy E 2 is defined as and the low mode energy E 1 is defined as where u = [u; v; w] is the velocity vector, ū is the cell average and the filtered velocity ũ is extracted from the filtered vector of DoFs for the solution variables: with û defined similarly and the filters σ , σ defined as a vector of size where k ∈ K are the polynomial modes and λ is determined iteratively as depicted in Figure 1.In other words, an Energy Ratio (ER) (Tantikul and Domaradzki 2010) is calculated based on the ratio of high mode to low mode energy; the high mode energy is the adaptively filtered solution Ũ subtracted by the strong filtered solution Û, which is defined as the zeroth-and first-order terms of the modal polynomial representation; the low mode energy is the strong filtered solution Ũ subtracted by the cell average Ū in order to maintain Galilean invariance.The adaptively filtered solution Ũ is the solution used for the numerical approximation and the expression in Equation ( 22) represents that the filter at minimum 1 for all modes and at maximum 0 for all modes n>1. Figure 1 shows that λ is incremented by 0.1 each time the criterion ER>UT is violated, with λ = 0 by default.The upper threshold, UT, was determined experimentally and is tabulated in Table 2.No grid/order independent normalisation factor was found, although it may be possible in a manner similar to that presented in Flad, Beck, and Munz (2016).This threshold represents the boundary between excessive filtering, which causes the simulation to crash, and insufficient dissipation, which leads the total kinetic energy to unphysically increase as seen in the unfiltered results presented in Figures 2 and 3.

Results
The method is tested on the nearly-incompressible Taylor-Green vortex, a popular choice for testing the dissipative properties of schemes in both well-resolved and under-resolved settings as well as, more generally, a simple model for turbulent flow and the cascade of energy from larger to smaller scale eddies.We examine the inviscid variant, a challenging problem due to the fact that the flow will always be under-resolved at late enough simulation times.The physical behaviour of the solution in the inviscid limit is still not fully understood; in particular, the existence of a finitetime singularity is debated (Majda and Bertozzi 2001).In addition, it is commonly believed that the total kinetic energy is preserved at all times for the inviscid Taylor-Green vortex; however, this is reliant on the assumption that there is an absence of singularities.Fehn et al. (2022) provides numerical evidence to the contrary, as well as a comprehensive review of theoretical and numerical analyses of the existence of singularities for three-dimensional incompressible Euler flows.
Overall, there are few published results on the inviscid Taylor-Green vortex to longer times, likely due to the difficulty of achieving robustness without being overly dissipative in the numerical method or oscillatory in the rate of kinetic energy dissipation, as well as the lack of DNS results.Shu et al. (2005) presents enstrophy and total kinetic energy results to t ≤ 6 using a Fourier collocation method with a sharp cutoff filter, a modal filter, and a WENO-FD method with resolutions up to 384 3 and shows that the method with the sharp cutoff filter is able to preserve total kinetic energy until the final time of the simulation.The kinetic energy dissipation rate is not shown.Colombo, Crivellini, and Nigro (2023) presents entropy and kinetic energy results to t ≤ 20 with an implicit, entropy conserving DG method in entropy variables and demonstrates that their scheme conserves entropy, and kinetic energy to a remarkable degree when an odd-order discretization is used, showing a decrease of < 15% for an 8 3 /P5 method and < 10% for a 32 3 /P3 method.The kinetic energy dissipation rate is also not shown.Gassner, Winters, and Kopriva (2016) is the most comprehensive, and presents kinetic energy, entropy, enstrophy, numerical viscosity, and kinetic energy dissipation rate results to t ≤ 14 using several robust split-form schemes with resolutions up to 256 3 .It is shown that the kinetic energy is not conserved for all times, and that the kinetic energy dissipation rate is slightly oscillatory despite the method being non-adaptive and peaks around 9s, although both the oscillations and overall amplitude is dampened when the resolution is increased.Chapelier, de la Llave Plata, and Lamballais (2016) presents kinetic energy dissipation rate results to t ≤ 20 using a modal DG method with overintegration, and we find similar results as shown in Figure 4 (Bustamante and Brachet 2012) presents results to t ≤ 5 with a resolution of 4096 3 , and using the Beale-Kato-Majde theorem with the analyticity-strip theorem, suggests that a finitetime singularity cannot be ruled out for the inviscid TGV.Fehn et al. (2022) presents results to t ≤ 20 using an incompressible method with resolutions up to 8192 3 , the highest resolution results to date, to the best of the authors' knowledge.These results suggest the existence of a singularity, as the kinetic energy is ostensibly grid-converged and reveal an onset of dissipation around t>5, while the maximum vorticity and enstrophy continue to increase at a consistent factor with increasing resolution.If there were no singularity, then the kinetic energy should be conserved until progressively later times as the resolution increases, and the vorticity supremum should approach a finite maximum.However, even at these resolutions, the possibility that the dissipation is purely numerical and not physical cannot be dismissed.No estimate for the required resolution to study the Taylor-Green vortex has been proposed, although (Cichowlas and Brachet 2005) suggests a resolution of ≈ 16-32 k is required to confirm the existence of a singularity for the Kida-Pelz flow, a variant of the TGV with additional symmetries.Therefore, it is likely that no numerical answer regarding the existence of a singularity in the inviscid TGV will be found in the near future.
The results presented in this paper are up to t ≤ 15 with resolutions up to ≈ 105 3  A hexahedral mesh is used to discretise the computational domain = [0, 2π] 3 with periodic boundaries and the initial condition which corresponds to a maximum Mach number of M ≈ 0.08 is given by the following profile of primitive variables: v(x, 0) = − cos(kx) sin(ky) cos(kz), ( 24) We plot the total kinetic energy k E , defined as: with respect to the time t; the kinetic energy dissipation rate, defined as the negative rate of change over time of the total kinetic energy; and the enstrophy, defined as: where ω is the vorticity vector.Several reference DNS results at varying Reynolds numbers from Brachet et al. (1983) and Fehn et al. (2022) are plotted in Figure 5.At sufficiently high Reynolds numbers, such as all the ones plotted, the kinetic energy dissipation rate peaks at some value around t ≈ 9, with the Re = 3000 and Re = 5000 results being nearly identical.Of particular note is the inviscid result from Fehn et al. (2022), in which the total kinetic energy is conserved until around t ≈ 5 at which point the dissipation begins and peaks in a similar manner to the viscous results where Re > 3000.Results with lower Reynolds numbers and/or spatial resolutions are omitted, but both factors in general lead to earlier onsets of dissipation and lower peaks.
The inviscid TGV case is run with an 8 3 and a 16 3 hexahedral mesh using the previously described modal DG method with polynomial orders 3, 4, and 5 used for the Taylor basis.In Figures 2 and 3, results for the unfiltered DG method are presented.As expected, almost all cases diverge at some point.The P3 method for both grids reach a peak between t ≈ 5 − 6 and smoothly decline to zero, while the higher order methods oscillate before crashing; this is likely due to the increased amount of smaller-scale structures being captured by the P4 and P5 methods.The 16 3 /P5 result is able to run fully.The onset of dissipation is earlier than in the (Fehn et al. 2022) result and there are some wiggles in the curve, but it is remarkably similar to the Re = 5000 result; further analysis is required for this result.
Figures 6 and 7 present the results using the previously described AF-DG method.The 16 3 /P5 method was not run with AF as it did not suffer from instabilities in the unfiltered version.The dissipation rate is noisy, a side effect of the adaptive filtering process.It can be seen that the filtering is engaged around t ≈ 5 in all cases, which is consistent with the development of the hypothetical singularity.The estimated peak of a naked-eye curve-fit is around t ≈ 6 for the P3 method on both grids; the P4 and P5 methods on the 8 3 grid peak around t ≈ 8, although it is difficult to tell due to the noisiness of the data; the 16 3 /P4 method after some oscillations peaks around t ≈ 9.This is promising as it suggests that at this low resolution, the overall behaviour of the flow is already somewhat approximated.As the number of DoFs is increased, the onset of dissipation is delayed, with the latest being t ≈ 4 for the 16 3 /P5 case, while the results of Fehn et al. (2022) show that the onset of dissipation should approach t ≈ 5 as the resolution increases.
The maximum value for λ, the filter strength, is plotted in Figure 8.Here again it can be seen that the filtering is engaged after t ≈ 5.For the P3 methods, the filter value is extremely low.For the P4 and P5 methods on the 8 3 grid, the amount of filtering done is relatively sparse, while for the 16 3 /P4 method, the filtering is much more frequent.It is noted that all methods do not diverge.In addition, the noisiness of the dissipation data does not show in the total kinetic energy plot.
When the mesh is increased to 32 3 , plotted in Figure 9, the P3 method runs in the unfiltered mode but has an incorrect dissipation curve, and the P5 method displays similar behaviour to the 16 3 /P5 results, showing the inefficiency of increasing mesh resolution rather than increasing the order of the numerical method.The P4 method crashes even with Figure 6.AF-DG TGV on an 16 3 hexahedral mesh, resolutions of ≈ 43 3 , 52 3 , 61 3 , and 8192 3 , respectively.AF.This warrants further investigation, but the most likely reason is the even-ordered method and the smaller difference between the grid resolution and the polynomial approximation degree, which decreases the effectiveness of filtering and is also mentioned in Flad, Beck, and Munz (2016).The filtering operation is only effective when there is a large range of scales within each cell.
Figure 10 plots the enstrophy of the AF-DG and the unfiltered 16 3 /P5, 32 3 /P3, and 32 3 /P5 methods; however, as discussed in Shu et al. (2005), the enstrophy is not physically meaningful when the flow is under-resolved.The results somewhat resemble those presented in Fehn et al. (2022), with the enstrophy increasing until t ≈ 7, at which point it tapers off.This is consistent with the hypothetical singularity: the enstrophy and maximum vorticity increase as the resolution increases until it can no longer be resolved and numerical dissipation takes over.In our implementation, the enstrophy of the 8 3 /P5 case which continues to rise as time increases suggest that the physics are incorrectly resolved for this case, and the 32 3 /P5 case does not display an increased amount of enstrophy as it should according to its resolution.Similarly, the  maximum vorticity at each timestep for the unfiltered 8 3 /P3 and 16 3 /P5 and the AF 8 3 /P3, 8 3 /P4, and 8 3 /P5 methods is plotted in Figure 11.For the converged cases, the maximum vorticity tapers off as the flow develops, whilst for the unfiltered 8 3 /P3, the maximum vorticity can be seen to diverge.These quantities are plotted semilogarithmically in Figures 12 and 13.For the enstrophy, the tapering-off is easier to see.For the maximum vorticity, the growth rate appears to be e 2 3 t until t ≈ 3.5, which is reported in Fehn et al. (2022) as well.
Figures 14 and 15 plot the Laplacian term D, the effective viscosity ν e , and a dimensionless function i , which are defined as in Aspden et al. (2008): ) and  where x is the cell length.i is representative of small-scale numerical dissipation, which is also related to the Kolmogorov length scale as η e = i x.At the 8 3 grid resolution, the P3 and P4 results display similar behaviour to that shown in Aspden et al. (2008) although the peak is earlier than expected, consistent with the other quantities shown in this work.The effective viscosity, as expected, decreases with increasing resolution as well as with the simulation time.The dimensionless function i peaks, then approaches a constant value.This is shown by the 8 3 /P4 result, although not by the P3 and P5 results, which suggest that the small-scale numerical dissipation is somewhat accurately captured by the P4 scheme by not the others.These quantities are representative of a scheme's numerical dissipation properties, and warrant further investigation.In particular, a filtering criterion based on these quantities could be superior to the presently used energy ratio threshold.The test cases were also run with N + 1 quadrature, which was enough to guarantee stability.The results are presented in Figures 16 and 17, and reveal a lower  dissipation for 8 3 /P4 and 32 3 /P4, the cause of which needs further investigation.However, the 8 3 /P3 and 16 3 /P4 results show promising dissipation behaviour, albeit with earlier onset due to the increased numerical dissipation compared to the 8192 3 result.The 32 3 /P4 case was cut short due to computing time restrictions.Overall, this shows that an increase of quadrature accuracy to N + 1 attains stability for our method without the use of filtering or any other stability mechanisms.However, despite the success in achieving stable solutions with over-integration at this resolution the development of the adaptive filtering method is still important because, as discussed in Moura et al. (2017), as the polynomial order increases, the instability is not caused by aliasing errors but rather the sharp dissipation characteristics of higher-order DG methods.

Conclusion
This paper presents an adaptive filtering modal DG method, which has been tested at varying resolutions and polynomial orders on the inviscid Taylor-Green vortex.The method was able to handle the aliasing instabilities that are often present in the high-order solution of under-resolved flows.The simulation of turbulent flow in under-resolved settings is still a developing topic.The results presented in this paper suggest that a no-model highorder DG method with relatively few degrees of freedom is able to resolve a significant amount of small-scale features.An adaptive dissipation method is developed which chooses an appropriate filter strength based on the ratio of high-mode to lowmode energy.The primary drawback is the lack of a filtering criterion which does not require tuning.This method allows for stability despite aliasing errors at low resolutions, but creates noisy dissipation behaviour.
The key uncertainty regarding simulations of the inviscid TGV is that of the physical solution itself.Despite the remarkable resolutions to which the problem has been solved, it is still uncertain whether or not a finite-time singularity develops.However, recent evidence in the literature has suggested that the answer is yes, as simulations with increasing resolutions return seemingly converged kinetic energy results while vorticity and enstrophy continuously increase.This would mean that the ability of numerical methods to solve the TGV until late times is solely due to numerical dissipation and does not reflect the physical solution of the problem.Therefore, it  is difficult to draw conclusions regarding the accuracy of the method based on inviscid TGV results; however, it provides a challenging test case against which the eddy-resolving, dissipation, and robustness characteristics of numerical methods can be tested.
A great deal of future work is required to improve the AF-DG method and advance our understanding singular or nearly-singular flows.Increasing the polynomial order allows for increasing the resolution of the solution more efficiently than through grid refinement.However, as previously discussed and reported in the literature, higher order DG schemes display sharper dissipation properties, which causes instability through a stronger bottleneck effect.Therefore, improved stability measures are required at higher orders.In addition, more analysis of high-order DG schemes, especially with filtering, is required.Determination of a grid-/order-independent criterion for the upper threshold of the energy ratio both improves the practicality of the method and may offer additional information about the underlying dynamics.Finally, studies should be extended to turbulence in viscous and wall-bounded flows.

Figure 7 .
Figure 7. Maximum filter strength for AF-DG TGV on 8 3 and 16 3 hexahedral meshes.The λ value is the number of polynomial modes filtered, as detailed in Equation (22).

Table 1 .
Williams, Shunn, and Jameson (2014)s for hexahedra, with the nodal quadrature points defined by the tensor product Gauss-Legendre points and the modal ones defined byWilliams, Shunn, and Jameson (2014).

Table 2 .
Upper threshold for energy ratio for various grid resolutions and polynomial orders, determined experimentally for the inviscid TGV case.