Resolution and Energy Dissipation Characteristics of Implicit LES and Explicit Filtering Models for Compressible Turbulence

Solving two-dimensional compressible turbulence problems up to a resolution of 16, 3842, this paper investigates the characteristics of two promising computational approaches: (i) an implicit or numerical large eddy simulation (ILES) framework using an upwind-biased fifth-order weighted essentially non-oscillatory (WENO) reconstruction algorithm equipped with several Riemann solvers, and (ii) a central sixth-order reconstruction framework combined with various linear and nonlinear explicit low-pass spatial filtering processes. Our primary aim is to quantify the dissipative behavior, resolution characteristics, shock capturing ability and computational expenditure for each approach utilizing a systematic analysis with respect to its modeling parameters or parameterizations. The relative advantages and disadvantages of both approaches are addressed for solving a stratified Kelvin-Helmholtz instability shear layer problem as well as a canonical Riemann problem with the interaction of four shocks. The comparisons are both qualitative and quantitative, using visualizations of the spatial structure of the flow and energy spectra, respectively. We observe that the central scheme, with relaxation filtering, offers a competitive approach to ILES and is much more computationally efficient than WENO-based schemes.


Introduction
The use of nonlinearly weighted and upwind-biased schemes (such as weighted essentially non-oscillatory (WENO) reconstructions) is considered a state-of-the-art numerical tool for the simulation of turbulent compressible flows in the presence of shock waves [1][2][3].These numerical schemes have an inherent numerical dissipation mechanism due to their asymmetric stencil reconstructions.They are also extensively used in the implicit large eddy simulation framework for turbulent flow computation and represent a promising approach to meeting the degree of accuracy required in high precision engineering requirements [4][5][6][7][8].On the other hand, non-dissipative central (symmetric) linear reconstruction schemes can be used in a relaxation filtering framework where an explicit low-pass spatial filtering procedure is utilized to remove high-frequency contents of the simulation [9][10][11].The main idea of this paper is to examine the relative advantages and disadvantages of these approaches for solving canonical turbulent compressible flow systems including a stratified shear layer instability problem to demonstrate the evolution of linear perturbations into a transition to nonlinear two-dimensional hydrodynamic turbulence.
The stratified Kelvin-Helmholtz instability (KHI) is a famous problem which manifests itself when there is a velocity difference at the interface between two fluids of different densities [12].It can commonly be observed through experimental observation, numerical simulation and it is also visible in many natural phenomena such as for example in situations with wind flow over bodies of water causing wave formation and in the planet Jupiter's atmosphere between atmospheric bands moving at different speeds [13].The study of this instability in a benchmark formulation reveals key information about the transition to turbulence for two fluids moving at different speeds.For these practical applications, it is common to choose a dual shear layer (DSL) setup to simulate the formation of KHI in a two-dimensional framework of the Euler equations for the validation of numerical methods [14].The evolution of the instability can be compared across simulations with different underlying numerical schemes qualitatively (through visual examination of variable contours) and through statistical tools such as the measurement of kinetic energy spectra, the total energy in the flow and the rate of dissipation of total energy.Our goal is to categorize the differences in the averaged kinetic energy spectral scaling capture of the different numerical methods used in this investigation and pronounce conclusions about their advantages and disadvantages.We also aim to examine field plots in the form of density contours to visually assess the extent of spatial filtering.Next we study another test case of a Riemann shock interaction problem given by [15] (hereby denoted RSI) which will be used primarily to assess the shock capturing ability of the investigated numerical methods.It must be noted that KHI type vortical structures are also observed in this test case [16].Through this investigation, we hope to connect the choice of solver to the choice of physics that needs to be captured for the numerical methods (and broader fluid mechanics) community.
The numerical methods utilized in this investigation are from two chief fields of thought (with regards to modeling approach).We use the implicit large eddy simulation methodology [5,[17][18][19][20][21][22] where weighted non-oscillatory (WENO) upwinding schemes [23][24][25] are coupled with certain Riemann solvers in an implicit large eddy simulation (or ILES) framework.The term implicit refers to the addition of numerical dissipation by the upwinding scheme rather than any explicit dissipation term [26].It is no surprise, therefore, that the characterization of the dissipative behavior of various WENO type schemes is an important area of interest in the computational fluid dynamics community [27][28][29].Although the primary motivation for the development of non-oscillatory upwinding schemes was for the simulation of hyperbolic conservation laws with sharp discontinuities in the solution field [30][31][32][33][34], ILES has emerged as a popular alternative to the more formal techniques for subgrid scale modeling in turbulent flows [35][36][37][38][39] since it is convenient to let the numerical dissipation inherent to these schemes substitute the subgrid scale dissipation that would otherwise require artificial modeling [40][41][42][43].We remark that controlling this numerical dissipation is not a trivial task, but the computational efficiency of this approach (as it does not require any explicit turbulence model) makes it attractive.
We shall also utilize the concept of relaxation filtering [44][45][46] (which employs low-pass spatial filters for removing high frequency content) and utilize a central scheme for the reconstruction of conserved quantities at cell interfaces.In this framework, the flow field variables are filtered after a certain number of timesteps to dissipate the amount of energy related to residual stresses and thus the dissipative effects of the unresolved scales on the resolved scales are modeled.We aim to utilize the dissipative behavior of this framework to prevent Gibbs oscillations which arise due to the use of central schemes in governing laws with dominating hyperbolic behavior.Relaxation filtering has been used with success for the large eddy simulation of homogeneous isotropic turbulence [47] and represents a computationally inexpensive approach to subgrid scale modeling [48].Various filters have also been utilized for the relaxation filtering methodology [49,50] and the choice of free filtering parameters is crucial for this method.In addition, the relaxation filtering approach has been observed to have the desirable quantity of being grid independent [51].
We note that the primary purpose of both these approaches is the addition of dissipation in coarse grained simulations so as to capture kinetic energy spectra scaling as well as the formation of shocks.We remind the reader that the addition of artificial dissipation is crucial for the accurate simulations of hyperbolic conservation laws.The aforementioned numerical methods are implemented using the message passing interface (MPI) approach for distributing a physical domain among multiple processes for concurrent computation [52].Through this, exceptional high fidelity simulations (HFS) are obtained for the DSL and RSI problems as benchmark 'true' numerical experiments for the purpose of determining the viability of a numerical method for coarse grained usage.On a thorough comparison of density contours, kinetic energy spectra capture and total energy measurement, conclusions are drawn about the suitability of numerical methods for the purpose of numerical simulation using the two-dimensional Euler equations.We reiterate that the main idea of this work is to provide the reader an exhaustive comparison of the well studied WENO-Riemann methods and the symmetric (central) reconstruction based explicit filtering approach through the aforementioned conclusions combined with a detailed discussion devoted to the mathematical development of these numerical approaches.The domain decomposition used in our MPI methodology is also reported briefly for the purpose of reproducibility.We aspire to provoke questions in our readers' mind about the performance of the numerical schemes discussed in representing accurate physics for two-dimensional flows given by the DSL and RSI test cases which are a simpler representation for many real world applications.For the convenience of the reader, we can draw up a list of points to summarize the major questions this investigation attempts to address: 1.
How does a central scheme augmented with a relaxation filter compare with the performance of the conventional upwind-biased approach to numerical dissipation in ILES? 2.
What reduction in computational expense can we expect through the utilization of this new explicit filtering framework?3.
Can we quantify the effect of the free modeling parameters that accompany low-pass spatial filters on the statistical and instantaneous features of our solution field? 4.
What parameters dominate the dissipative nature of a WENO-ILES scheme for our test cases?
The rest of our work will be devoted to the development of the mathematics and a documentation of results to address these important questions.

Governing Equations
As mentioned in the previous section, we shall utilize the two-dimensional Euler equations in their conservation form as our underlying partial differential equations.These can be expressed as [53,54] ∂q ∂t where where Here ρ, P, u and v are the density, pressure, horizontal and vertical components of velocity respectively.Also, e and H denote the total energy and total enthalpy per unit mass with γ being the ratio of specific heats.Before proceeding further, we develop the eigensystem of the system of equations which will then be used in devising hyperbolic conservation laws.The convective flux Jacobian matrices for our system of equations are and where . A similarity transform of the above flux Jacobian matrices can be shown as where Λ and Ψ are the diagonal matrices comprising of the real eigenvalues of A and B respectively; L and S represent 4 × 4 matrices whose columns are the eigenvectors of this eigendecomposition.They can be expressed as [55] where a is the speed of sound and can be calculated through the relation a 2 = γP/ρ.The parameter β is given by 1/(2a 2 ).We must reiterate that this is one of several possible eigendecompositions possible for this system and there may be slight variations in results for certain algorithms according to one's choice [56].The solvers we evaluate in this paper are independent of the choice of eigenvector matrices.

Finite Volume Framework
The semi-discrete form of the governing equations can be written as with q i,j being the cell-averaged vector of dependant variables, F i±1/2,j representing the fluxes at the right and left cell boundaries and G i,j±1/2 representing the fluxes at the top and bottom boundaries of the cell.We use the method of lines to represent our system of PDE's as an ODE through time (so as to implement a Runge-Kutta scheme for time integration) with the right hand side of the above equation representing the combined effect of all the spatial derivatives in the governing equations.The above ODE representation can be advanced to obtain the solution field at a future time step n + 1 given the current time step n by a total variation diminishing third-order Runge Kutta scheme [57] as follows q (1) where a time step ∆t is prescribed through a CFL criterion as In our work, we have defined the parameter η = 0.5 for all simulations.

WENO Reconstructions
The left and right states of the cell boundaries are reconstructed using the weighted essentially non-oscillatory (WENO) approach which was first introduced in [58].The order of accuracy of these reconstructions depends on the length of the stencil chosen and affect the solution (through a different dissipative behavior).In this work, we shall focus on the 5th order accurate WENO scheme (WENO-5) which utilizes a five point stencil.The reader is directed to [14] for a thorough examination of the effect of WENO stencil sizes on the Kelvin-Helmholtz instability test case.In what follows (and for the rest of this document), we shall introduce stencil expressions only for the x direction of our domain.We note, for the purpose of clarity, that it is the conserved variables which are being reconstructed in this formulation.We utilize a modified implementation of the WENO-5 reconstruction [59] which can be given as where the nonlinear weights are defined by with the smoothness indicators defined as 14) 15) In this study, we compute the smoothness indicators for each conserved variables.The optimal linear weighting coefficients are d 0 = 1/10, d 1 = 3/5 and d 2 = 3/10 in Equation ( 11) and d 0 = 3/10, d 1 = 3/5 and d 2 = 1/10 in Equation (12).The presence of any discontinuity thus leads to an adaptation in the weights for order reduction and increased dissipation in the stencil corresponding to the discontinuity.This expression for the calculation of the nonlinear weights of the smoothness indicator for the WENO-5 interfacial flux reconstruction will be referred to as WENO-JS.We also investigate the effect of choosing a different approach for the calculation of the nonlinear weights (hereby denoted as WENO-Z) which are given by [60] Before proceeding, we must observe that there are many variants of the WENO reconstructions described above (for example in [61][62][63][64]) and this investigation is by no means exhaustive.The different variants of the WENO reconstructions have all been developed as a response to demands of either different orders of accuracy [65,66] or reduced computational expense.The interested reader is directed to [33] for an excellent discussion on the effect of the order of reconstruction on a classical two-dimensional configuration utilizing the Euler equations.
We remark, for clarity, that the left (L) and right (R) states correspond to reconstructions considering the possibility of advection from both directions.The value of is set to 10 −6 for the WENO-JS approach and 10 −20 for the WENO-Z approach in order to prevent any errors by division of zero.We note that the power of the denominator in α k is usually set to p = 2 in most utilizations of the WENO-5 approach.In this investigation, we aim to perform a sensitivity analysis for p to quantify its effect on our solutions.A mathematical analysis of the role of the smoothness indicator and the parameter can be found in [67].

Riemann Solvers
Once the left and right states have been constructed using WENO-5, we may now use Riemann solvers to calculate the fluxes at these boundaries.We utilize four celebrated Riemann solvers for our numerical simulations which are as follows.We emphasize that the selection of these solvers was motivated by the desire to span a wide range of dissipative behaviors as reported in literature [16].

Rusanov Scheme
The Rusanov solver [68] utilizes information from the maximum local wave propagation speed to give us the follow expression for the flux where F R and F L are the flux components using the right and left constructed states respectively (i.e., F R = F(q R i+1/2 ) and F L = F(q L i+1/2 )).The local wave propagation speed c i+1/2 , is given by the maximum absolute eigenvalue of the Jacobian matrix of F between the cells i and i + 1, i.e., where r(A) represents the spectral radius of matrix A. For our case, the spectral radius may simply be determined as r(A) = max(|u|, |u − a|, |u + a|).Thus our final expression for the wave propagation speed in a cell can be written as

Roe Scheme
The Godunov theorem states the the exact values of the face fluxes can be computed by the following relation if the Jacobian matrix is constant [69] where |Λ| is the diagonal matrix obtained after a similarity transformation of the Jacobian matrix A. For the Euler equations however, A is a function of the conserved variables (i.e., A = A(q)).
To account for this, the Roe Riemann solver [70,71] estimates the interfacial fluxes in an approximate Godunov manner where the superscript tilde represents the Roe density weighted average between the left and right states.Our modified eigensystem can be constructed by using the following averaged values of the conserved variables instead where the left and right states of the un-averaged conserved variables are available from the WENO-5 reconstruction described earlier.We note that the Roe averaged speed of sound is now a function of these averaged variables.We also utilize Harten's [72] entropy fix approach for expansion shocks by replacing Roe averaged eigenvalues with where ã is the Roe averaged speed of sound and λi are the eigenvalues of the Roe averaged Jacobian matrix Ã.Here is a small positive number and is typically chosen as = 0.1 in our simulations.

HLL Scheme
Harten et al. [73] assumed that the lower and upper bounds on the characteristic speeds can be used in the solution of the Riemann problem.These upper and lower bounds were given by [74][75][76] where S L and S R are the lower and upper bounds on the left and right state characteristics speeds.
The HLL solver thus takes the following form

AUSM Scheme
The advection upstream splitting method (or AUSM) was proposed to provide an alternative approach to low-diffusion flux-splitting methods [77,78].This approach is based on recognizing that the inviscid flux consists of a distinct convective flux and a pressure flux.A cell interface advection Mach number is defined to determine convective flow quantities.Our investigation uses the AUSM approach for its lower dissipative behavior through the following expression of the interface numerical flux where in which the directional convective Mach number (M = u/a) is given as The following split formula is used for the pressure We note that there are many variants of the classic AUSM solver [79][80][81][82][83][84], which are used to obtain more accurate and robust results in all-speed regimes.

Central Reconstruction Schemes
This approach relies on a traditional finite volume formulation for calculating the density, energy and velocities at the cell interfaces following which the cell face fluxes are calculated as a function of these reconstructed variables.The following 6-point stencil symmetric non-dissipative scheme is used for face reconstruction of the conserved quantity [85] q i+1/2 = a(q i+1 + q i ) + b(q i+2 + q i−1 ) + c(q i+3 + q i−2 ) (34) where the stencil coefficients are given by Once the relevant face quantities are determined from the nodal values, the fluxes may be calculated for use in Equation (7).The considerable simplicity of this approach is apparent.This simple central scheme may be combined with a classic relaxation filtering method as shown below to remove high wavenumber content and prevent floating point overflow.This can be shown to be an explicit large eddy simulation approach for conservation laws [44].

Relaxation Filtering
A low-pass spatial filter sweep is carried out using a seven point stencil (as given in [86]), wherein it is denoted as a standard relaxation filtering procedure if carried out a the end of time integration or as a modified relaxation filtering procedure when carried out at each substep of the TVDRK3 algorithm.The motivation of low-pass spatial filtering is to remove high frequency content from the solution field at every sweep and thus prevent any accumulation of noise (which might eventually lead to Gibbs oscillations).For the purpose of filtering, we use the following expression where and σ = [0, 1] is a parameter that controls filter dissipation strength with σ = 0 being completely non-dissipative and σ = 1 being most dissipative.The most dissipative σ = 1 case corresponds to a Gaussian-type filter with a full attenuation at the smallest grid cut-off scale.The choice of applying this filter sweep at each substep of the time integration or at its conclusion would appear to have ramifications on its success in preventing oscillations in the solutions for the conserved variable.This shall be examined in detail in the results section.As expected, the use of this explicit LES approach for the discrete hyperbolic equations leads to a significant reduction in computational expense.However, it may be prone to stability issues when used for strongly shocked flows (particularly for lower values of the free parameter controlling dissipation).It must be noted here that there is no restriction on the choice of filtering strategies when it comes to the relaxation filtering approach, any low-pass spatial filtering process (of course with a similarly shaped transfer function) would achieve similar results [87].The transfer function of the filter chosen here is plotted in Figure 1 and can be represented as where ω is the modified wavenumber.A value of σ = 1 gives us the 6th order binomial smoothing filter which is derived in [88].Using a Taylor series expansion procedure in Equation (36), it is easy to show that the leading truncation error term in this linear filter becomes sixth order accurate where h is considered our grid discretization size and the Roman numeral superscript on q denotes the 6th order differentiation.We note that Equation ( 36) is used in each direction sequentially (i.e., first in y direction, the direction of domain decomposition, to eliminate any extra MPI data exchange cost due to the filtering processes).

Shock Filtering
The shock filtering approach devised by Bogey et al. [89] relies on an adaptive filtering approach to remove noise from the solution field.The adaptive quantum of dissipation is calculated from the local solution to modify the free parameter (σ) as mentioned in the previous subsection.The basic kernel of the shock-capturing filtering approach is where In this work, we choose n = 1 to give us the following simplified version of the above equation with c 0 = 0.25 and c 1 = −c 0 .The filtering parameter σ s is calculated adaptively through extracting the high-wavenumber components of the pressure field using a second-order filtering process given by pi = (−p i+1 The magnitude of the high-pass filtered pressure squared term is then calculated as We can now devise a shock sensor through the ratio r given by where s = 10 −16 is needed to avoid numerical divergence.We note that it is possible to determine a shock detection through the local dilatation rather than the pressure (which is particularly useful for discerning turbulent fluctuations and shocks), but we limit the investigation to the pressure magnitude version of the shock sensor to align ourselves with the WENO-Riemann solver architecture which have primarily been developed for the purpose of shock capturing.Our filtering parameter may now be updated dynamically by For r i ≤ r th , the filtering magnitude σ s = 0. Values of r i > r th give us nonzero values of σ s with σ s → 1 for r i → +∞.The threshold parameter r th is thus our free parameter for the addition of dissipation and has been set between values of 10 −4 and 10 −7 .Our cell face value for the filtering parameter σ s can be approximated as We remark here that this nonlinear explicit filtering scheme is expected to be more dissipative than the standard (or modified) central relaxation filtering approach in the previous subsection since it employs a considerably lower order of reconstruction (2nd order).

Results
This section is devoted to a quantification of the performance of the various numerical approaches outlined above.We first define the initial conditions of the two test cases before outlining our MPI methodology followed by which the results of the different numerical solvers are explained in the form of density contours, kinetic energy spectra and total energy plots.

Problem Definitions
As mentioned previously, we utilize the DSL problem and the RSI test case for the evaluation of our numerical methods.The DSL problem consists of the simulation of an initially perturbed compressible shear layer in two dimensions.This test problem characterizes the solvers' ability to evolve a sine wave perturbation into two-dimensional turbulence.We expect to see the instability trigger small-scale vortical structures at the sharp density interface initially which eventually transition (through nonlinear interactions) to a completely turbulent field.We must note here that the turbulent field obtained after sufficient progress in the simulation are due to growth of the shear layer instabilities.Figure 2a describes the initial conditions for the DSL problem.Our computational domain is a square of unit side length with periodic boundary conditions in all directions.The initial conditions are specified as follows One can observe that the vertical component of the velocity is perturbed using a single mode sine wave (n = 2, L = 1.0) and the amplitude of the perturbation is set at λ = 0.01.The DSL numerical experiments are solved to a final dimensionless time of t = 5.The second test case we have selected for this investigation is that of a shocked flow (i.e., the RSI test case) given by configuration 3 in [15] as shown in Figure 2b.Details about the configuration in question may also be found in [16].This problem is designed with different initial conditions specified in each quadrant of our computational domain.These initial conditions may be varied to obtain a variety of combinations of shocks and contact discontinuities.The final time we consider for the RSI problem is t = 0.5 and the transmissive boundary conditions are used in both directions (e.g., q b−k/2 = q b+k/2 , where k = 1, 2, and 3 and the subscript b refers to the boundary edge).In both problems, the ratio of specific heats is considered to be constant γ = 7/5.The R 2 simulation domain for all experiments is set (x, y) ∈ [−0.5, 0.5] × [−0.5, 0.5].
Both test cases are evaluated using three coarse grid resolutions given by N 2 = 256 2 , 1024 2 and 4096 2 points with a N 2 = 16, 384 2 solution obtained using the Roe Riemann solver as a reference high fidelity simulation (HFS).Density contours and kinetic energy spectra scaling are provided for the purpose of evaluation of the numerical schemes in the DSL test case.Classical k −3 scaling from KBL theory is also provided for the purpose of comparison [90][91][92].Only density contours are provided for the RSI test case (for visually ascertaining shock capturing ability).

MPI Methodology
For the purpose of parallelization, the OpenMPI [52] standard library is used to speed up our computation.Using the familiar principles of distributed multiprocessing, our physical domain is split up into several smaller computational domains which are then distributed to different computational units (or processes).Each process solves the governing equations on its own slice of the physical domain before communicating boundary data to its neighbors.This domain decomposition and data transfer can be observed in Figure 3.One can observe that there is only one direction in which data exchange is needed, since the physical domain is distributed into even slices with the same number of points in the x direction.It must be noted that in case of the DSL problem, data must also be exchanged between the top and bottom computational slices cut from the physical domain due its periodic boundary conditions.Data exchange is programmed to occur at the end of each time integration substep as well as during the explicit filtering process.This is because the calculation of derivatives and the filtering process requires data values from beyond the computational domain allocated to the process.This extra data belonging to the neighboring domains is stored in 'ghost' arrays for each computational slice.The ghost arrays must thus be updated concurrently with the solution through send and receive commands.The process of sending and receiving data from these ghost arrays is done using non-blocking send MPI_Isend and non-blocking receive MPI_Irecv commands.These non-blocking communication privileges imply that the function calls return before the communication is completed.This is done to ensure that there is no deadlock.A deadlock implies a concurrent use of synchronous send or receive buffers (which may only be used by one thread).The use of the aforementioned non-blocking commands alleviates this issue but also requires the use of the MPI_wait command to ensure all communication is complete before the code executes any further on each thread.The interested reader is referred to [93] for a up to date tutorial on the MPI standard.A computational performance assessment of our code is detailed in Figure 4 where both strong and weak scalings are described.In a weak scaling test, the problem size is proportionally scaled with the number of processes, so that the computational work per core remains essentially constant.Ideally, the relative runtime (the run time scaled by that on one process) should remain close to unity as the number of processes and the problem size are increased [94].The weak scaling test demonstrates how the parallel overhead scales as a function of NP and provides some measure of how the communication complexity of an algorithm degrades as the process count increases.In a strong scaling test, a fixed problem size is considered and a varying number of processes are used to solve the problem.The parallel method should exhibit run times inversely proportional to the number of processes used, or alternatively, the speedup should vary linearly with NP.The strong scaling test will also show limiting behavior for a fixed problem size, since for a large enough NP, there will not be sufficient computational work available on each process to adequately overlap parallel communication overhead.Speedup (S) is defined as the ratio of execution time of the serial algorithm to execution time of the parallel algorithm with NP processes (i.e., S = t 1 /t NP ).
One can observe a marked increase in relative runtime (and also a corresponding decrease in efficiency) as the order of the processes (NP) becomes closer to the order of the number of grid points in each direction.This is because the relative cost of data exchange becomes higher.From a computational perspective, a future improvement of the current implementation of our code would be to decompose the physical domain through a resolution dependent heuristic for optimal efficiency.We note, however, that for the purpose of this investigation, the speed up offered enabled us to obtain our reference HFS data (i.e., approximately 268 million degrees of freedom) in approximately 400 hours of simulation time (on 432 processes) for the DSL case and 50 h (on 288 processes) for the RSI case.

Dual Shear Layer (DSL) Problem
In this subsection, we perform an exhaustive investigation into the performance of the different numerical techniques described previously on the DSL test case.Figure 5 shows two density contours (at t = 5) for the Roe and HLL Riemann solvers at the highest resolution (i.e., corresponding to N 2 = 16, 384 2 ).The high degrees of freedom leads to the capture of several small scale structures in both solvers.We remark here that a considerable loss in correlation can be observed between the two solver solutions manifesting itself in a phase difference.This is due to the different amount of dissipation imparted to the system by each Riemann solver -the different perturbations cause a different evolution of the flow field (a hallmark of deterministic chaos).We remark here that the research in numerical scheme development should ideally be aligned to the goal of improved predictions of both statistical quantities such as kinetic energy spectra and of phase.Figure 5 also describes the total energy evolution of the system.It must be noted that an ideal evolution of the Euler equations would conserve energy for all time but the dissipation due to our numerical methods eventually reduce the total kinetic energy of the system.Once the instability transitions to a turbulent field, one can see a loss of phase between the HLL and Roe solvers from the total energy evolution plot.The total kinetic energy of the solution field is defined as the spatial average of the instantaneous kinetic energy at all points and our averaging operator is given by where A is the area of the physical domain and E(x, y, t) is the instantaneous kinetic energy per unit mass at a particular point in the solution field (i.e., E = 1/2(u 2 + v 2 )).Indeed, we compute the integral using the cell centered values of velocity components.Figure 5 also details the difference between the HLL and the Roe solvers in terms of final time kinetic energy spectra.It is clear that the flow fields (for these HFS) have shown the k −3 scaling from KBL theory in the inertial range.The existence of this scaling has also been observed in other benchmark two-dimensional configurations (e.g., a single-mode Richtmyer-Meshkov instability problem [95]) using the same underlying governing equations.From the comparison of the HLL and Roe Riemann solvers in terms of spectral scaling, one can visually represent the difference in the dissipative nature of both solvers.We observe, as expected, that both approaches represent the critical inertial and integral length scales with the same degree of accuracy.The kinetic energy spectra that we have represented here can be calculated using the following expression in wavenumber space [96] Ê where û(k, t) is the Fourier transform of the velocity vector in the wavenumber space.Equation ( 54) can be also rewritten in terms of velocity components where we compute velocity components û(k, t) and v(k, t) using a fast Fourier transform algorithm [97].Finally, the spectra can be calculated by angle averaging in the following manner where k = |k| = k 2 x + k 2 y .Figure 6 shows the time evolution of the solution field in the form of three snapshots of the density contours at t = 1, 3 and 5 for the Roe Riemann solver.We have also included the transition to turbulence for our reference HFS.These snapshots are also provided for the coarser grids and show a considerable reduction in the small scale details with increasing coarseness.For the coarsest resolution of N 2 = 256 2 one can see that there are very few small scale details and only large scale structures are preserved.The absence of smaller structures with increasing coarsening also gives evidence of spatial filtering being performed implicitly by the ILES approach.We can also observe the presence of acoustic waves in the central belt of high density in the earlier stages of the instability.Similar trends are seen in Figure 7 which details the coarse grid simulations for the HLL Riemann solver.The slightly higher dissipation of this approach is apparent with a reduced number of smaller structures for a similar resolution when compared to the Roe approach (a visual examination of the N 2 = 1024 2 case confirms this).This can also be seen in Figure 8 where the Rusanov solver is used for the construction of interfacial fluxes.Both the Rusanov and HLL Riemann solvers are considered amongst the more dissipative approaches for solving a system of hyperbolic equations.The Rusanov approach happens to be the most dissipative solver we have investigated here.This becomes clear when the kinetic energy spectra for each solver are compared (later).Figure 9 shows the performance of the AUSM solver (which is the least dissipative Riemann solver in this investigation) and it can be seen that it triggers a turbulent phase (due to a reduced dissipation of numerical noise) even in the coarsest of runs.Spectral scaling confirms this lack of dissipation.(a (a The DSL problem was also simulated using the central scheme for constructing interfacial fluxes and subsequent low-pass spatial filtering (dubbed the relaxation filtering approach).This approach produced accurate solutions for the DSL problem without the use of any upwinding as seen in Figure 10 for a free filtering parameter value of σ = 0.5.It can also be observed that small scale structure capture is similar to the performance of the Roe solver.This is an important observation-the proposed relaxation filtering approach now allows us to control the dissipation in an active manner through its filtering parameter.This can be observed through density contours for a lower value of σ = 0.25 as shown in Figure 11 where much more smaller structures are captured and for the highest value of σ = 1.0 where larger scales are preserved at the expense of the smallest ones (in Figure 12).This potentially allows for experimentation with different classes of low-pass spatial filters (with different transfer functions) in order to obtain the desired amount of scale capture.We also note the considerable computational ease of implementation for this framework and reduced computational expense in comparison to the relatively involved nature of the WENO reconstructions and the Riemann solvers.We make a note of caution though that the choice of the filtering parameter and the associated filter also play a role in the stability of the numerical scheme [98].It is entirely possible that a very low quantity of dissipation may lead to the emergence of nonphysical Gibbs oscillations and possible floating point overflow (i.e., σ = 0 refers to a perfectly central scheme without any numerical dissipation).In this regard, the Riemann solvers used with the WENO-5 reconstruction are decidedly superior with their conditional upwinding using nonlinear weights in their stencils.One possible approach to improve the stability of this approach is to implement the low-pass spatial filtering sweep at the end of each time integration substep.In theory, this would improve the stability of the relaxation filtering approach.We find that the relaxation filtering done at multiple substeps performs just as well as the regular relaxation filtering approach and it shall also be seen that this approach does not represent a considerable computational overhead in comparison to the regular relaxation filtering approach.The density contours of this modified implementation of the relaxation filtering approach are shown in Figure 13.To the authors' experience, the relaxation filtering approach with the filtering strength of σ = 0.5 yields accurate results on both fine and coarse grids in practice, although it doesn't fulfil full attenuation at the cut-off grid scale (i.e., see Figure 1).
Figure 10.DSL density contours with time obtained through the standard relaxation filtering approach (using σ = 0.5) with coarse grid resolutions of (a (a (a Figure 13.DSL density contours with time obtained through the modified relaxation filtering approach (using σ = 1.0) with coarse grid resolutions of The adaptive implementation of the relaxation filtering approach (called the shock filtering approach here) was also implemented to solve the DSL problem as shown in Figure 14.Here, as mentioned previously, the control parameter for dissipation is the shock threshold parameter r th which controls the sensitivity of the shock sensor.At this high value of r th = 10 −5 , it is seen that the shock sensor provides the least amount of dissipation and results in turbulence being triggered even for the coarsest runs.In fact, on closer examination, spurious numerical oscillations can also be seen in the density contours.It is also observed that the larger vortices in the coarse grained simulations are less well defined.This fact will be confirmed later in the kinetic energy spectra comparisons where we shall also observe a deviation from the expected inertial range scaling.On reducing the magnitude of r th we increase the dissipation of the shock filtering scheme as shown in Figures 15 and 16.It is seen that for coarser runs, the increased dissipation manages to damp the initial perturbation to the flow field and prevent any sort of transition to turbulence (as seen in the case with r th = 10 −7 .
(a (a (a Figure 16.DSL density contours with time obtained through the shock filtering approach (using r th = 10 −7 ) with coarse grid resolutions of Figure 17 shows the performance of the aforementioned numerical schemes quantified in a statistical sense through kinetic energy spectra scaling at different resolutions for the four different WENO-Riemann combinations.Plots showing the total energy in the flow are also reported.It is generally seen that the flow transitions to turbulence somewhere around t = 1.For the purpose of comparison, we utilize HFS statistical data from the Roe Riemann solver (with N 2 = 16, 384 2 ).We also include the expected theoretical scaling given by KBL theory for kinetic energy cascades in two-dimensional turbulence (i.e., E(k) ∝ k −3 ).This figure details the dissipative behavior of the ILES techniques for three different coarse grid simulations.It can be seen that the integral length scales are captured well by each WENO-Riemann combination whereas the difference in dissipative behavior becomes apparent at the higher wavenumbers in the cascade.Our previous observations about the dissipation of small scale structures by the different Riemann solvers in the ILES approach are confirmed where it can be seen that the AUSM solver is seen to be the least dissipative in comparison to the other approaches.The Rusanov solver is seen to be the most dissipative.The Roe and HLL approach are generally seen to be in between the Rusanov and AUSM approaches.This is also confirmed through the total energy content plots which show the Rusanov and HLL solvers reducing the total energy content of the flow at the same rate until turbulence ensues.We detail the dissipative behavior of the standard relaxation filtering approach (i.e., relaxation filtering at the end of time integration) in Figure 18.The effect of the filtering parameter σ is clearly seen with higher values of σ proving more dissipative.We note that using a different family of filters would perhaps lead to a greater control over the dissipative characteristics of this approach.The interested reader is directed to [99] for a comparison of the effects of different filters on large eddy simulation type approaches.In general, a very good approximation of the HFS kinetic energy spectra is obtained by this method.The integral length scales are also captured relatively well.The total energy content plots also confirm the dissipative behavior of the explicit filter with larger values of σ generally reducing the total energy content at a faster place.A curious observation here is that the relaxation filtering approach destabilizes the flow at a comparatively similar point in time for all coarse resolutions.This could imply a resolution independent dissipation phenomenon.Figure 19 shows the performance of the shock relaxation filter for the DSL problem.It can immediately be noticed that the integral length scales are captured incorrectly at the coarsest resolutions by this approach.The highest value of the shock threshold (r th = 10 −4 ) fails to return a solution (due to floating point overflow) for the N 2 = 4096 2 resolution but albeit with limited success in terms of scaling capture for the other two resolutions.This points to a resolution dependent dissipation mechanism which would imply a need to vary the shock sensor threshold manually.Increasing the resolution of our simulations leads to a better performance at the integral length scales by this approach but a considerable proneness to aliasing error is seen at the cut-off wavenumbers (particularly for higher values of the shock threshold parameter).A comparatively larger amount of inertial range dissipation is also seen in the form of a downward shifting of the coarse grained kinetic energy spectra.This can be attributed to the lower order filtering procedure when compared to the higher order stencil used in the standard relaxation filtering approach.The total energy content plots also outline the same trends with far more rapid fall off in the total energy of the solution field in comparison to the ILES and standard relaxation filtering method.We clarify that the r th = 10 −4 simulation eventually led to a floating point overflow for the N 2 = 4096 2 case (as seen in the final time kinetic energy spectra plot).We tentatively conclude that the proposed relaxation filtering mechanism performs in a better manner than the adaptive shock capturing approach for this investigation.
Figure 20 shows the effect of the WENO-5 flux reconstruction parameter p on both versions of our nonlinear weight calculator (i.e., WENO-JS and WENO-Z) for the coarse resolution of N 2 = 1024 2 .The kinetic energy spectra and total energy plots shown here are for the Roe Riemann solver.An inspection of both these quantities proves inconclusive about the effect of p when it comes to the trends in dissipation.On examining Figure 21, we can see once again that the choice of p does not affect the kinetic energy spectra at the final time of the simulation.However, the total energy content plots show that lower values of p preserve the total energy of the flow before it becomes turbulent.Figure 22 shows the effect of the choice of the flux reconstruction (i.e., WENO-JS or WENO-Z) for the ILES approach.It is seen that the choice of the nonlinear weight calculation assumes more importance when utilizing a more dissipative solver (such as the HLL or Rusanov Riemann solver) where the WENO-Z reconstruction proves to be less dissipative than the WENO-JS approach.On the contrary, for the Roe and AUSM solvers (which are less dissipative), both methods perform in the same manner.A final examination of the DSL problem is made with respect to the implementation of the relaxation filtering approach.Kinetic energy spectra from the normal implementation (where the relaxation filtering sweep is carried out at the end of each time integration) and the modified implementation (with filtering at the end of each timestep) are compared in Figure 23 for four different values of the filtering parameter σ at the resolution of N 2 = 1024 2 .We observe that both approaches add the same amount of dissipation for integral length scale capture of features and show slight differences in dissipative behavior at the highest wavenumbers.We note that spurious oscillations can be arrested with much more certainty in the modified implementation of the relaxation filtering kernel.The relative dissipative behavior of both implementations remains relatively constant over the scale of σ values.

Riemann Shock Interaction (RSI) Problem
The RSI test case is also studied for the shock capturing ability of the previously described numerical schemes.Density contours at t = 0.5 are displayed for each of the different numerical approaches for the purpose of determining their relative differences in dissipation and shock capturing ability.Figure 24 shows the HFS density contours obtained using the four Riemann solvers combined with WENO flux reconstructions.At this fine resolution, it is seen that all four WENO-Riemann combinations are equally able to capture the smallest relevant length scales in the flow.The shock evolution of this configuration is generally expected to be symmetric about the center line of the shock plume and the same is noticed in the Roe, HLL and Rusanov solvers.The AUSM Riemann solver, however, shows a slight skew in the formation of the plume which may be associated with low dissipation inherent to the scheme.It is possible that this low dissipation coupled with the Kelvin-Helmholtz vortices being formed during the evolution of the shock lead to the observed skewing.Figure 25 shows the results of the coarse grid resolution runs for this test case using the 4 Riemann solvers, it is apparent from the visual examination of the density contours shown here that the AUSM solver remains the least dissipative of all the Riemann approaches studied here and shows the skewing of its shock even at the relatively coarse grid run of N 2 = 4096 2 .A successive reduction in the presence of small scale features can also be seen with increasing coarseness in our simulations.The coarsest runs also display non-physical striations (commonly known as shock oscillations) in the density contours which are due insufficient dissipation at that level of grid size.These striations are most pronounced for the Roe and AUSM solvers as one might expect.These shock oscillations are commonly visible in the contour plots when their lower frequency modes remain undamped.A characteristics based reconstruction procedure (i.e., interpolating the characteristic variables) would reduce these oscillations, albeit with added expense.Figure 26 shows the performance of the standard relaxation filtering employed at the end of the time integration procedure (using a filter parameter of σ = 0.5) for three coarse resolutions as well as performance of shock relaxation filter for three different r th values.It can be seen that the standard relaxation filtering approach leads to a good result for the density field and may be considered comparable in quality with the shock capturing WENO-Riemann schemes.The shock filtering approach is seen to be adept at capturing the most sharp gradients in the field but higher higher values of r th lead to a loss in the symmetry of the shock plume with the evidence of a skewing (for example in subfigure (g)) even without small scale capture.The highest value of the shock capturing threshold parameter used r th = 10 −5 displays a considerable skewing of the shock plume with much less small scale capture as compared to the standard relaxation filtering kernel with σ = 0.5.It must be observed here that the striations observed in the ILES runs do not appear to be visible using this methodology.This is due to the low-pass spatial filtering characteristics of the explicit filtering kernel.The higher frequency components of these oscillations however may be present in the solution field as they are usually present to a certain extent in each computational realization even on very high fidelity simulations.Figures 27 and 28 compare the WENO-JS and WENO-Z approaches to cell face flux reconstruction for the RSI test case for two different values of p. From the density contours displayed at a coarse resolution of N 2 = 1024 2 and 4096 2 , it is difficult to conclude about the dissipative behavior of either reconstruction mechanisms for this particular case.However, the faint appearance of striations for the WENO-JS case with p = 2 seems to imply that this particular configuration of the reconstruction mechanism is less dissipative in comparison to the others.Finally, Figure 29 shows a one to one comparison of the instantaneous density contours for the relaxation filtering approach in its standard (i.e., low-pass spatial filtering at the end of the time integration) and modified (i.e., low-pass spatial filtering at each substep of the time integration) approach with the filtering strength of σ = 0.5.A visual examination does not show any large deviation in the dissipative behavior of both implementations of the relaxation filtering approach with the absence of any major skewing of the shock plume.However an examination of the N 2 = 1024 2 instance does indicate a more dissipative behavior of the modified approach.We note however, the difference in the dissipative behavior of either of these implementations is minimal (as shown in the double shear layer test case through the examination of kinetic energy spectra).The added number of filter sweeps acts as a preventive methodology for shock discontinuities leading to Gibbs phenomena.Table 1 shows a compilation of computational expense for the RSI test case for the three coarsest resolutions each using a different number of processes.We have used 12 process for N 2 = 256 2 , 48 processes for N 2 = 1024 2 and 192 processes for N 2 = 4096 2 .A broad message from this table can be obtained from the final column which tabulates CPU times for N 2 = 4096 2 case where the relaxation and shock filtering approaches are seen to be around 4 times faster (on average) than the WENO-Riemann solver ILES method.This adds value to the conclusions stated above, wherein the standard relaxation filtering approach is seen to perform as well as the ILES approach but at a much lower computational expense.We must also note that the modified standard relaxation approach does not add a significant amount of overhead to the standard relaxation filtering methodology.We must clarify that the tabulated times at the lower resolutions are not indicative due to the shorter duration of simulations and the unequal load distribution on the shared HPC cluster.Indeed, the Rusanov solver requires the least computational work load due to its simplicity.It does not involve any conditional statements or characteristic transformations.One possible reason for this inhomogeneous performance of compute clusters could be inhomogeneous hardware equipment (i.e., uneven CPU or memory equipment of the nodes).Another reason could be software related interactions.We report here our assessments only for a particular simulation at a certain time.Running the same executable on a different occasion with a different set of compute nodes would yield results in a slightly different wall-clock time.Therefore, their assessments will be valid only in a statistical sense.However, we would like to highlight that the central scheme, with relaxation filtering, offers a competitive approach to ILES and is much more computationally efficient than WENO-based schemes.

Conclusions
This study presents an in depth analysis of both implicit large eddy simulation (ILES) methods and explicit filtering approaches for the purpose of simulating compressible hydrodynamic turbulence by solving the two-dimensional Euler equations.We start by outlining the mathematical formulation of our governing equations followed by which we introduce the two major computational approaches we have used: (i) upwind-biased ILES formulations and (ii) central non-dissipative reconstructions equipped with explicit filtering procedures.To test the fidelity of the numerical methods being examined here, we use the dual shear layer (DSL) test case as well as a Riemann shock interaction (RSI) problem to evaluate their kinetic energy spectra resolving and shock capturing capabilities.Both test cases are simple representations for the physics of many real world phenomena and are well suited for benchmarking the advantages and disadvantages of our chosen numerical schemes.Evaluations of the different numerical schemes are done through the comparison of statistical data such as the kinetic energy spectra and total energy as well as through qualitative analyses of the instantaneous flow structures.A final comparison of the computational expense of the different numerical approaches is also reported.
In order to examine the general nature of the upwind-biased based numerical schemes, we utilize four different WENO-Riemann solver combinations (i.e., the Roe, HLL, AUSM and Rusanov solvers) with slightly different characteristics in terms of dissipative behavior and fine scale structure capture ability.The upwinding required for using these Riemann solvers in their implicit LES formulation is obtained through nonlinear biased interfacial flux constructions using the fifth order accurate WENO approach.Two versions of the smoothness indicators (the WENO-JS and WENO-Z formulations) are used for the reconstruction and the parameter p in these smoothness indicators is also varied for a sensitivity analysis.These ILES solvers (with different flux reconstruction and p configurations) are used to simulate coarse grained solution fields with N 2 = 256 2 , 1024 2 , and 4096 2 resolutions for the DSL and RSI problems.These coarse grained runs are compared to high fidelity simulation data obtained through the simulation of the test cases at N 2 = 16, 384 2 using the Roe Riemann solver.The theoretical scaling from KBL theory is also provided for the purpose of comparison.Predictably, it is observed that the Roe and AUSM solvers are less dissipative as compared to the HLL and Rusanov solvers and that all implementations of the Riemann solvers in the ILES framework are successful in capturing the integral length scales (along with a fair length of the inertial range in the turbulent cascade).Simulations which utilize different values of p are seen to remain more or less invariant in terms of kinetic energy spectra scaling.It is noted, however, that the choice of the smoothness indicator (i.e., whether WENO-JS or WENO-Z) has important ramifications on the dissipative nature of the more dissipative HLL and Rusanov Riemann solvers for which the WENO-JS scheme is seen to be more dissipative.The low-dissipation of the AUSM Riemann solver also causes a skewing of the shock plume generated in the RSI test case when compared to the other Roe, Rusanov and HLL approaches.In the RSI case, nonphysical striations (i.e., shock oscillations) are also observed at the coarsest simulation of N 2 = 256 2 for the WENO-Riemann solver combinations.These striations are most pronounced, once again, for the less dissipative AUSM approach.In our future studies, we plan to examine the results of the WENO interpolations done using the characteristics variables and in terms of their effects on these oscillations.
The DSL and RSI test cases are also simulated using the relaxation filtering approach where a low-pass spatial filtering sweep is applied to the solution field.This sweep may be applied at the end of each timestep (which we call the standard relaxation filter) or at each substep of the time integration process (denoted the modified relaxation filter).The dissipative nature of the filter is controlled by the free parameter σ which varies between 0 and 1 with the former implying no filtering and the latter denoting highest dissipation.We can remark here that the σ = 1 case corresponds to a standard Gaussian shaped filter.A shock capturing version of the aforementioned filter is also used where the filtering strengths are calculated adaptively through a heuristic based on the local solution field.The model parameter which controls dissipation is a shock sensor r th which is varied from 10 −7 to 10 −4 with lower values being more dissipative.It is seen that the standard relaxation filtering approach performs in a manner similar to the ILES approach with excellent integral length scale capture as well as improved inertial range scale capture.It is also shown that a significant decrease in computational expense is also obtained.The standard and modified relaxation filters are also shown to be more or less similar to each other in dissipative behavior.The modified relaxation filter is implemented and tested primarily due to the fact that the added number of filter sweeps at each timestep would be less susceptible to Gibbs oscillations.In addition to the standard and modified relaxation filtering approach, the shock filtering methodology is tested on the DSL and RSI test cases for different values of its shock threshold parameter where it is observed that integral length scales are captured poorly at the coarsest resolutions.It is seen that this filtering technique is prone to over-dissipating in the inertial range (which shows itself as a downward shift of the averaged kinetic energy cascade).A pile up of kinetic energy at the cutoff wavenumber is also seen for the shock threshold values tested in this work.Although a proper choice of r th may result in an accurate prediction, we conclude that the determination of the optimum value for the threshold parameter r th is not that trivial and depends on the resolution as well.The explicit filters were also applied to the RSI test case where it was seen that the low frequency shock oscillations visible in the N 2 = 256 2 test case were dissipated effectively by the choice of our free parameters.However, it was noted that the shock filtering approach caused a distinct skewing of the shock plume at higher values of r th .Once again, the performance of the standard and modified relaxation filters was comparable to the ILES approach.
We conclude by revisiting the questions posed in the introduction and may summarize our findings as follows: 1.
It can be stated with a certain degree of confidence that the central schemes implemented with relaxation filters perform in a comparable manner to the WENO-Riemann solver combinations tested here and that the standard and modified relaxation filters represent a viable alternative to the ILES framework for the test cases investigated here.Not only are they seen to capture small scale features and their nonlinear interactions accurately, they also appear to perform well when it comes to capturing shocks.However, WENO based ILES schemes are attractive, because they do not require any ad-hoc tuning or filtering parameters.

2.
We note a significant reduction in computational complexity due to the relative simplicity of the numerical algorithm for the relaxation filtering framework.In addition, the computational expense of this framework is seen to be much lower in comparison to the different WENO based ILES solvers.We find that the relaxation filtering (or explicit filtering) approach is approximately four times faster than the WENO reconstruction based ILES approaches.

3.
The free modeling parameters associated with the proposed numerical framework perform as expected with a changing transfer function shape corresponding to a modification of the dissipative behavior of the numerical method.This also allows us to actively control the dissipation in the numerical method through model free parameters σ and r th for the standard and shock relaxation filtering approach respectively.

4.
We have performed a thorough investigation of the WENO-Riemann solver based ILES schemes where it is verified that the AUSM Riemann solver is the least dissipative whereas the Rusanov solver is the most.The flux reconstruction scheme does have an effect on the dissipative behavior of the ILES methodology when it comes to the expressions for the calculation of the smoothness indicators (i.e., WENO-JS and WENO-Z), but for the dissipative Riemann solvers (Rusanov and HLL) only.The WENO-JS scheme is seen to be more dissipative than the WENO-Z scheme for these two solvers and more or less identical for the less dissipative AUSM and Roe solvers.It is seen that p does not have a major effect on the dissipative effect of the scheme for this investigation.However, less dissipative results are obtained for lower p values.
Our closing comments are made on the fact that the standard or modified relaxation filtering framework could possibly represent a viable tool for the efficient computation of hyperbolic conservation law problems.The main conclusion is that the central scheme, with relaxation filtering, offers a competitive approach to ILES and is much more computationally efficient than WENO-based schemes.We aim to perform further investigations for other benchmark flows to lend more strength to this conclusion.

Figure 1 .
Figure 1.Transfer function for low-pass spatial filter used in standard and modified relaxation filtering.Notice reduction of area under curve with increasing values of σ implying more dissipation.

Figure 2 .
Figure 2. Representative initial conditions for the (a) dual shear layer (DSL) and (b) Riemann Shock Interaction (RSI) test cases.Note that we use periodic boundary conditions on all sides for the DSL problem while using open boundary conditions in all directions for the RSI test case.

Figure 3 .
Figure 3.A schematic of our domain decomposition.Data is exchanged in the y direction prior to each time integration substep and during the explicit filtering process.

Figure 4 .
Figure 4. Message passing interface (MPI) performance indicators.(a) Weak scaling keeps the problem-per-node size constant while increasing the number of nodes.(b) Strong scaling spread the same size problem across more nodes.

Figure 5 .
Figure 5. High fidelity (N 2 = 16, 384 2 ) density contours at t = 5 for the Roe (a) and HLL (b) Riemann solvers.Accordingly, the evolution of the total kinetic energy content (c) as well as the kinetic energy spectral scaling (d) of the flow are shown for both cases.Note the phase difference in the solutions once it transitions to turbulence.

Figure 6 .
Figure 6.DSL density contours with time obtained through the Roe solver.High fidelity data with N 2 = 16, 384 2 (a-c) and coarse grid simulations at N 2 = 4096 2 (d-f), N 2 = 1024 2 (g-i) and N 2 = 256 2 (j-l).Note absence of smaller structures in the coarse grid runs implying spatial filtering being performed implicitly.

Figure 17 .
Figure 17.Kinetic energy spectra (a,c,e) for the DSL problem obtained using four different weighted essentially non-oscillatory (WENO)-Riemann solver combinations at time t = 5. High fidelity kinetic energy spectra from a Roe simulation at N 2 = 16, 384 2 data points included as reference.Total energy of the flow (b,d,f) also plotted from t = 0 to t = 1.

Figure 18 .
Figure 18.Kinetic energy spectra (a,c,e) for the DSL problem obtained using the standard relaxation filtering approach at time t = 5. High fidelity kinetic energy spectra from a Roe simulation at N 2 = 16, 384 2 data points included as reference.Total energy of the flow (b,d,f) also plotted from t = 0 to t = 1.

Figure 20 .Figure 21 .
Figure 20.Effect of p in the smoothness indicator on dissipative behavior for the DSL problem.Sensitivity is checked with both the WENO-JS and WENO-Z implementations for the Roe solver at N 2 = 1024 2 with kinetic energy spectra (a,b) and total energy plots (c,d).

Figure 24 .
Figure 24.A comparison of the effect of different WENO-Riemann solver combinations on density contours at t = 0.5 for the RSI test case at a resolution of N 2 = 16, 384 2 .(a) Roe; (b) HLL; (c) AUSM; (d) Rusanov.

Figure 27 .Figure 28 .
Figure 27.A comparison of the (c,d) WENO-JS and (a,b) WENO-Z scheme (for different p values) through effect on density contours for the RSI test case at a resolution of N 2 = 1024 2 for the Roe solver.

(a) N 2 =Figure 29 .
Figure 29.A comparison of the different implementations of the central scheme with the relaxation filtering procedure with σ = 0.5 (i.e., standard (a,c,e) or modified (b,d,f)) for the RSI test case.

Table 1 .
CPU times for different solvers for the RSI test case.12 processes were used for N 2 = 256 2 , 48 processes were used for N 2 = 1024 2 and 192 processes were used for N 2 = 4096 2 .