Comparative statistics of selected subgrid-scale models in large eddy simulations of decaying, supersonic MHD turbulence

TL;DR

This study compares traditional and recent subgrid-scale models in large eddy simulations of decaying, supersonic MHD turbulence, finding the nonlinear model with a small explicit filter best captures higher order statistical features.

Contribution

It evaluates the effectiveness of EV, SS, and NL SGS models in compressible MHD turbulence, highlighting the nonlinear model's superior performance with small explicit filters.

Findings

Nonlinear SGS model improves higher order statistics.

Explicit filtering influences SGS model impact.

No statistical influence from models when using grid-scale quantities.

Abstract

Large eddy simulations (LES) are a powerful tool in understanding processes that are inaccessible by direct simulations due to their complexity, for example, in the highly turbulent regime. However, their accuracy and success depends on a proper subgrid-scale (SGS) model that accounts for the unresolved scales in the simulation. We evaluate the applicability of two traditional SGS models, namely the eddy-viscosity (EV) and the scale-similarity (SS) model, and one recently proposed nonlinear (NL) SGS model in the realm of compressible MHD turbulence. Using 209 simulations of decaying, supersonic (initial sonic Mach number of ~3) MHD turbulence with a shock-capturing scheme and varying resolution, SGS model and filter, we analyze the ensemble statistics of kinetic and magnetic energy spectra and structure functions. Furthermore, we compare the temporal evolution of lower and higher order statistical moments of the spatial distributions of kinetic and magnetic energy, vorticity, current density, and dilatation magnitudes. We find no statistical influence on the evolution of the flow by any model if grid-scale quantities are used to calculate SGS contributions. In addition, the SS models, which employ an explicit filter, have no impact in general. On the contrary, both EV and NL models change the statistics if an explicit filter is used. For example, they slightly increase the dissipation on the smallest scales. We demonstrate that the nonlinear model improves higher order statistics already with a small explicit filter, i.e. a three-point stencil. The results of e.g. the structure functions or the skewness and kurtosis of the current density distribution are closer to the ones obtained from simulations at higher resolution. We conclude that the nonlinear model with a small explicit filter is suitable for application in more complex scenarios when higher order statistics are important.

Figures

Tables

Table 1. Table 1: Filter weights for a discrete representation of the box filter based on an optimal filter approach by minimizing the residual between the analytic and discrete filter over wavenumbers below the filter width. The filter width Δ ¯ ¯ Δ \overline{\Delta} is given in terms of grid-spacing Δ x subscript Δ 𝑥 \Delta_{x} .

Identifier	filter width	filter weights
	$\overline{\Delta }$	$w_i$	$w_{i\pm 1}$	$w_{i\pm 2}$
GS	$1{\Delta }_x$	$1$	$0$	$0$
F3	$2.711{\Delta }_x$	$0.4015$	$0.29925$	$0$
F5	$4.7498{\Delta }_x$	$0.20238$	$0.22208$	$0.17673$

Table 2. Table 2: Effective initial kinetic and magnetic Reynolds numbers in the simulations depending on resolution. The numbers are estimated according to  Zhou  et al. ( 2014 ) .

Resolution	$\mathrm{Re}$	$\mathrm{Rm}$
${128}^3$	$𝒪\left(600\right)$	$𝒪\left(450\right)$
${256}^3$	$𝒪\left(1100\right)$	$𝒪\left(900\right)$
${512}^3$	$𝒪\left(1700\right)$	$𝒪\left(1400\right)$

Table 3. Table 3: Computational efficiency of the different SGS models relative to the no-model run at identical resolution (ILES-128). For reference, the efficiency of the no-model run at the next higher resolution (ILES-256) is also shown. The numbers represent the mean values (at the fluid level) over three test runs of 300 cycles each. Each run used the same executable on a single machine employing 8 MPI-processes (no threading). This corresponds to a 64 3 superscript 64 3 64^{3} grid per process, as suggested by the Enzo documentation.

	GS	F3	F5
NL	1.175(5)	1.468(5)	2.411(7)
EV	1.133(5)	1.431(6)	2.368(8)
SS		1.949(8)	4.804(16)
ILES-256	9.242(28)111 Please note that this number only includes the time per cycle. The total computational costs are increased by another factor of $\gtrsim 2$ due to the decreased timestep at higher resolution.

Equations

$$\frac{\partial \overline{\rho }}{\partial t}+\nabla \cdot \left(\overline{\rho }\tilde{\mathbf{u}}\right)=0,$$

$$\begin{array}{r}\frac{\partial \overline{\rho }\tilde{\mathbf{u}}}{\partial t}+\nabla \cdot \left(\overline{\rho }\tilde{\mathbf{u}}\otimes \tilde{\mathbf{u}}-\overline{\mathbf{B}}\otimes \overline{\mathbf{B}}\right)+\nabla \left(\overline{P}+\frac{{\overline{B}}^2}{2}\right)\\ =-\nabla \cdot \tau ,\end{array}$$

$$\frac{\partial \overline{\mathbf{B}}}{\partial t}-\nabla \times \left(\tilde{\mathbf{u}}\times \overline{\mathbf{B}}\right)=\nabla \times \mathcal{E},$$

$${\tau }_{ij}={\tau }_{ij}^{\mathrm{u}}-{\tau }_{ij}^{\mathrm{b}}+\left(\overline{B^2}-{\overline{B}}^2\right)\frac{{\delta }_{ij}}{2},$$

$${\tau }_{ij}^{\mathrm{u}}\equiv \overline{\rho }\left(\widetilde{u_i{u}_j}-{\tilde{u}}_i{\tilde{u}}_j\right)\text{and}{\tau }_{ij}^{\mathrm{b}}\equiv \left(\overline{B_i{B}_j}-{\overline{B}}_i{\overline{B}}_j\right).$$

$$\mathcal{E}=\overline{\mathbf{u}\times \mathbf{B}}-\tilde{\mathbf{u}}\times \overline{\mathbf{B}}.$$

$$\overline{E}=\underset{(\mathrm{resolved})}{\underbrace{\frac{1}{2}\overline{\rho }{\tilde{u}}^2+\frac{1}{2}{\overline{B}}^2}}+\underset{=E_{\mathrm{sgs}}^{\mathrm{u}}+E_{\mathrm{sgs}}^{\mathrm{b}}\equiv E_{\mathrm{sgs}}(\mathrm{unresolved})}{\underbrace{\frac{1}{2}\overline{\rho }\left(\widetilde{u^2}-{\tilde{u}}^2\right)+\frac{1}{2}\left(\overline{B^2}-{\overline{B}}^2\right)}}$$

$${\hat{\tau }}_{ij}^{\mathrm{u}}$$

$$\hat{\mathcal{E}}$$

$${\nu }^{\mathrm{u}}$$

$${\eta }_t$$

$${\hat{E}}_{\mathrm{sgs}}^{\mathrm{u},{\mathcal{S}}^{\ast }}$$

$${\hat{E}}_{\mathrm{sgs}}^{\mathrm{b},\mathcal{M}}$$

$${\hat{\tau }}_{ij}^{\mathrm{u}}$$

$${\hat{\tau }}_{ij}^{\mathrm{b}}$$

$$\hat{\mathcal{E}}$$

$${\hat{\tau }}_{ij}^{\mathrm{u}}$$

$${\hat{\tau }}_{ij}^{\mathrm{b}}$$

$$\hat{\mathcal{E}}$$

$${\overline{\square }}_{i,j,k}=\sum _{l=-N^{\ast }}^{N^{\ast }}\sum _{m=-N^{\ast }}^{N^{\ast }}\sum _{n=-N^{\ast }}^{N^{\ast }}{w}_l{w}_m{w}_n{\square }_{i+l,j+m,k+n},$$

$$\mathrm{Re}=-⟨{\left|\nabla \times \mathbf{u}\right|}^2⟩L^2{\left(\frac{\mathrm{d}{\mathbf{u}}^2}{\mathrm{d}t}\right)}^{-1}.$$

$$\mathrm{Rm}=-⟨{\left|\nabla \times \mathbf{B}\right|}^2⟩L^2{\left(\frac{\mathrm{d}{\mathbf{B}}^2}{\mathrm{d}t}\right)}^{-1}.$$

$${\tan }^2\theta =\frac{\sum k_{\perp }^{2}E(k)}{\sum k_x^{2}E(k)},$$

$$\mathrm{skew}x=\frac{⟨{\left(x-⟨x⟩\right)}^3⟩}{{\sigma }^3\left(x\right)},$$

$$\mathrm{kurt}x=\frac{⟨{\left(x-⟨x⟩\right)}^4⟩}{{\sigma }^4\left(x\right)}-3.$$

$$S_{\parallel }^p(l)=⟨{\left|(\mathbf{u}(\mathbf{x}+\mathbf{l})-\mathbf{u}(\mathbf{x}))\cdot \mathbf{l}/l\right|}^p⟩$$

$${\zeta }_p=\frac{p}{9}+2\left(1-{\left(\frac{2}{3}\right)}^{p/3}\right).$$

Full text

I Introduction

Magnetohydrodynamic (MHD) turbulence is observed in many different processes and on many different scales, for example, in astrophysics, in magnetized accretion disks (Balbus and Hawley, 1998), stellar winds (Goldstein et al., 1995), galaxies and galaxy mergers Rodenbeck, Kai and Schleicher, Dominik R. G. (2016), or more generally in processes like magnetic reconnection (Oishi et al., 2015) and the turbulent amplification of magnetic fields (Schober et al., 2015). Moreover, experiments on Earth also aim at a better understanding of e.g. flow-driven MHD instabilities (Cooper et al., 2014). However, the full multi-dimensional dynamics are only rarely accessible in these observations and experiments. For this reason, simulations are frequently used as a third, complementary approach or to support the design of experiments (Tzeferacos et al., 2015).

Simulations of pure turbulence are nowadays possible at very high resolution Kaneda et al. (2003); Yeung et al. (2015); Federrath et al. (2016) and properly capture a lot of physical processes. However, there are still many extreme regimes, for example in astrophysics, where turbulence is thought to play an important role but which are inaccessible to direct simulations with realistic parameters. This situation is also not going to change in the near future despite the ever increasing availability and performance of large computing clusters, and the advances in numerical methods. In these cases, large eddy simulations (LES) have been employed successfully in the past, however, mostly in the realm of (incompressible) hydrodynamics, see e.g. (Sagaut, 2006) for a general introduction and (Schmidt, 2015) for an astrophysics related review. In LES only the largest eddies, which correspond to motions on large and intermediate scales, are simulated directly. The smallest scales, which are either not represented or unresolved in these simulations, are reintroduced by the means of a subgrid-scale (SGS) model. In other words, LES try to incorporate effects on the large scale flow that stem from the small scales or from interaction between large and small scales. Ideally, the overall quality of the simulation improves with respect to the physical processes that are taken into account.

Formally, this is equivalent to applying a low-pass filter to the ideal compressible MHD equations resulting in expressions of the form (Vlaykov et al., 2016)

[TABLE]

for a static, homogeneous and isotropic filter under periodic boundary conditions. The filtered primary quantities, i.e. the density $\overline{\rho}$, velocity $\widetilde{\bm{u}}$, magnetic field $\overline{\bm{B}}$ (incorporating $1/\sqrt{4\pi}$) and pressure $\overline{P}$ are considered resolved in LES. Normal filtering is denoted by an overbar $\overline{\Box}$ whereas mass-weighted, Favre (Favre, 1983) filtering is denoted by a tilde $\widetilde{\Box}=\overline{\rho\Box}/\overline{\rho}$.

Assuming an isothermal equation of state ($\overline{P}\propto\overline{\rho}$), all interactions between resolved and unresolved scales and among unresolved scales themselves are captured by the two new terms in the equations. The turbulent stress tensor is given by

[TABLE]

and can be decomposed into the turbulent (or SGS) magnetic pressure (last term), SGS Reynolds stress $\tau_{ij}^{\mathrm{u}}$ and SGS Maxwell stress $\tau_{ij}^{\mathrm{b}}$ with

[TABLE]

The second new term is the turbulent electromotive force (EMF):

[TABLE]

Moreover, the total filtered energy is given by

[TABLE]

where, by virtue of the identity $\tau_{ii}^{\Box}\equiv 2E_{\mathrm{sgs}}^{\Box}$ (with Einstein summation), the SGS energies are given by the traces of the corresponding stress tensor.

All these SGS terms are unclosed because the mixed terms, i.e. $\widetilde{u_{i}u_{j}}$, $\overline{B_{i}B_{j}}$ and $\overline{\bm{u}\times\bm{B}}$, are not explicitly accessible in an LES and require modeling. This modeling is the main challenge for a successful LES.

SGS models have been subject of research in hydrodynamics for many decades, e.g. in the incompressible (Sagaut, 2006), compressible (Garnier et al., 2009) and astrophysical (Schmidt, 2015) regime. However, work in the realm of MHD and in particular compressible MHD is scarce, see (Chernyshov et al., 2014) and (Miesch et al., 2015) for recent reviews. Directly linked to this work are the MHD simulations of (decaying) turbulent boxes in 3D (Miki and Menon, 2008), in 2D (Theobald et al., 1994) and in the incompressible case (Müller and Carati, 2002; Yokoi, 2013). However, all these groups use different numerical schemes, e.g. finite-differences or (pseudo-) spectral methods, while we employ a shock-capturing finite-volume scheme. Usually, these shock-capturing methods are thought to provide an implicit SGS model as shocks are captured by locally increasing the effective numerical dissipation with the help of e.g. slope limiting. Since this procedure is part of the overall method, these simulations are also referred to as implicit LES (ILES) Grinstein et al. (2007).

In the presented work, we compare the a posteriori behavior of several SGS model including a nonlinear model that explicitly captures compressibility Vlaykov et al. (2016). Their performance was previously evaluated a priori (Grete et al., 2015, 2016). Here, we analyze a set of simulations of decaying homogeneous and isotropic turbulence with respect to a set of statistical quantities including energy spectra, structure functions and statistical moments of the primary fields. More details on the numerics and the implementation are given in the following section II where we also introduce the different models tested and the setup of the particular test case of decaying compressible MHD turbulence. Afterwards, in section III we present the results of different statistics such as energy spectra, evolution of statistical moments and structure function. Then we discuss these results with respect to previous work in IV and conclude in the last section V.

II Method

II.1 Subgrid-scale models

In our a priori analysis (Grete et al., 2016) we tested three different model families, eddy-dissipation, scale-similarity, and nonlinear models, with different normalizations. All models were tested against the expressions (5) and (6) where the filtered nonlinear term is not split into additional components. A split allows to separate different types of interactions between scales: only among resolved scales, between resolved and unresolved scales, and only among unresolved scales. Thus, the previous a priori results evaluated the performance of the models to recover all kind of interactions simultaneously. This is also what we implicitly expect of the models in this work. For each model family we identified the best model with corresponding coefficients. These models are now tested a posteriori and described in the following.

The eddy-viscosity

($\mathrm{EV}$) model has the longest tradition with roots going back even further than its formulation for LES by Smagorinsky (Smagorinsky, 1963). While originally developed for the kinetic SGS stress tensor in hydrodynamics, the general idea has been transferred to MHD (Yokoi, 2013; Miki and Menon, 2008), where the EMF closure is usually referred to as anomalous or eddy-resistivity. The names of these functional models stem from their primary feature: purely dissipative behavior analogous to e.g. molecular viscosity and resistivity. We use the following formulations

[TABLE]

with the resolved traceless kinetic rate-of-strain tensor $\mathcal{\widetilde{S}}^{*}_{ij}=1/2\left(\widetilde{u}_{i,j}+\widetilde{u}_{j,i}\right)-1/3\delta_{ij}\widetilde{u}_{k,k}$ and current density $\overline{\mathbf{J}}=\nabla\times\overline{\bm{B}}$. The strengths of the eddy-viscosity $\nu^{\mathrm{u}}$ and eddy-resistivity $\eta_{t}$ are given by

[TABLE]

They are locally scaled by the SGS energies derived from realizability conditions of the SGS stresses resulting in (Vreman et al., 1994; Grete et al., 2015)

[TABLE]

with $\overline{\mathcal{M}}_{ij}=1/2\left(\overline{B}_{i,j}+\overline{B}_{j,i}\right)$ being the resolved magnetic rate-of-strain tensor. In agreement with the a priori analysis (Grete et al., 2016) the coefficients are set to $C_{1}=C_{2}=0.05$ and $C_{3}=C_{4}=0.04$. Please note that the coefficient indices are different from the ones in the referenced paper due to a reduced selection of models. Furthermore, the SGS Maxwell stress is neglected as the eddy-diffusivity extension in compressible MHD was found to not match the reference data (correlations $<0.1$) in the a priori analysis. Therefore, the effects of the SGS Maxwell stress are modeled implicitly by the numerical scheme.

The scale-similarity

($\SS$) model originates in experimental hydrodynamics where it was observed that scale-to-scale energy transfer is self-similar up to intermittency (Bardina et al., 1980). Formally, this additional scale separation can be expressed by an additional (test) filter \mathop{\vbox{\halign{#\cr\kern 3.0pt\cr\bracelu\leaders\vrule height=0.39998pt,depth=0.0pt\hfill\braceru$\crcr\kern 3.0pt\nointerlineskip\cr$\hfil\displaystyle{\Box}\hfil\crcr}}}\limits whose filter width is larger than the original one. The model is given by

[TABLE]

and mass-weighted filtering also applies to the test filter where velocity components are involved. Again, we choose the coefficients according to the a priori analysis: $C_{5}=0.67$, $C_{6}=0.9$ and $C_{7}=0.89$. The model allows for energy transfer down- and up-scale and, as a structural closure, aims at reproducing closely the properties of the SGS terms and not just their effects on the large scales.

The nonlinear

($\mathrm{NL}$) model is another structural model and exhibited the highest correlations with reference data in a priori tests (Grete et al., 2015, 2016). It can be derived from Taylor expansion of the inverse filter kernel (Yeo, 1987; Vlaykov et al., 2016) and requires no further assumptions about the underlying flow features. We employ the primary compressible extension resulting in the following model:

[TABLE]

As previously shown, this model does not require a calibration coefficient a priori Grete et al. (2016) and the prefactor $1/12$ originates from the second moment of the Gaussian or box filter.

II.2 Implementation and explicit filtering

We implemented the different models in the open-source, community-developed magnetohydrodynamic code Enzo (Bryan et al., 2014). The new terms are handled by operator-splitting within the MUSCL-Hancock framework and evaluated with centered finite-differences. They are advanced in time together with the other fluxes by the existing second order Runge-Kutta scheme.

Furthermore, we implemented a flexible filtering framework that supports different stencils and weights in real space. In order to determine the individual weights, we construct discrete, explicit filters by minimizing the residual between analytic and discrete filter yielding so called optimal filters (Vasilyev et al., 1998). We optimize for wavenumbers below the filter width (Sagaut and Grohens, 1999). The resulting weights and filter widths for a symmetric one-dimensional 3-point (F3), and 5-point (F5) stencil of the box kernel are listed in table 1.

We construct the corresponding multi-dimensional filters of $N$-point one-dimensional stencils by

[TABLE]

with discrete filter weights $w_{i}$, $N^{*}=(N-1)/2$, and indices referring to discrete spatial locations. This translates to the sequential application of the one-dimensional filters and results in large 3-d filter stencils ($N^{3})$, e.g. 27 points for the F3 filter and 125 points for the F5 filter. However, this construction is more accurate (Sagaut and Grohens, 1999) than the alternative approach of simultaneous application. Finally, we also use the trivial grid filter (GS). In that case, the quantities are used as they are computed in the original simulation, which corresponds to a natural filter by the discretization itself.

II.3 Simulations

All our simulations are run with Enzo (Bryan et al., 2014) using the HLL Riemann solver within the MUSCL-Hancock framework with second order Runge-Kutta time integration. Moreover, we close the equations of ideal MHD with a quasi-isothermal equation of state, i.e. we employ an ideal equation of state with a ratio of specific heats $\kappa=1.001$.

In order to get proper initial conditions for freely decaying, compressible MHD turbulence, we first start from uniform initial conditions $\rho=1$, $\bm{u}=\bm{0}$ and $\bm{B}=\left(0.6325,0,0\right)^{\mathrm{T}}$ (in code units) corresponding to an initial plasma beta (ratio of thermal to magnetic pressure) of $\beta_{\mathrm{p}}=5$. These uniform initial conditions are then driven in a cubic box with length $L=1$ and resolution $512^{3}$ by a stochastic forcing field that evolves in space and time (Schmidt et al., 2009). The forcing field has a parabolic profile between wavenumber $1<k<3$ and is centered at $k=2$. Moreover, the overall forcing amplitude is set to $V=3$ and distributed between $1/3$ compressive and $2/3$ solenoidal components. The forcing leads to statistically isotropic, homogeneous turbulence with root mean square sonic Mach number of $\mathrm{M_{s}}\approx 3$ after two turnover times $T=L/(2V)$ and we follow its evolution for a total of $20T$.

Assuming that two different snapshots are statistically independent from each other after one turnover time, we have an ensemble of 19 different realizations (at $t=\{2,3,\dots,20\}T$) of isotropic, homogeneous turbulence, which are statistically indistinguishable. We take states from this ensemble as initial conditions for freely decaying turbulence. This later allows us to analyze ensemble statistics to better capture the intermittent nature of turbulence.

For each realization we run the following 11 simulations with different configurations – varying the resolution, models (or lack thereof) and explicit filter, namely:

• •

3 implicit large eddy simulations (ILES). Recall that in these simulations there is no explicit model ($\tau={\bm{\mathcal{E}}}=0$), at resolutions of $128^{3}$, $256^{3}$ and $512^{3}$. They are referred to as ILES-128, ILES-256, and ILES-512, respectively.

• •

3 LES with the eddy-viscosity model at a resolution of $128^{3}$: one with no explicit filter (EV-GS); one with a filter with a three-point stencil (EV-F3); and one with a five-point stencil (EV-F5).

• •

2 LES with the scale-similarity model at $128^{3}$ employing three- (SS-F3) and five-point (SS-F5) explicit filtering (because a grid-scale scale-similarity model does not exist).

• •

3 LES with the nonlinear model again at $128^{3}$ with all three different filters NL-GS, NL-F3, and NL-F5.

The highest resolution $512^{3}$ ILES simulations are later used as reference runs. Comparing the results between the different ILES enables us to evaluate the pure influence of reduced resolution (and thus reduced dynamics) on the evolution of the decay. For example, the Reynolds number that can be achieved in a simulation depends on the resolution. Given that we use a shock-capturing scheme to solve the ideal MHD equations, i.e. viscosity and resistivity are of numerical nature rather than explicit, the kinetic and magnetic Reynolds numbers are not readily accessible. In ILES of decaying, homogeneous isotropic turbulence the effective kinetic Reynolds number can be estimated as Zhou et al. (2014)

[TABLE]

Extending this concept to obtain an estimate of the magnetic Reynolds number yields

[TABLE]

The resulting initial effective Reynolds numbers for the individual resolutions are listed in table 2. Thus, the effective magnetic Prandtl number of $\mathcal{O}\left(1\right)$. Finally, we can draw conclusions on the influence of the individual models on the evolution from a comparison of results between ILES and the different LES at identical resolution (here $128^{3}$)

Each simulation follows the decay for two turnover times. We capture snapshots every $0.05T$ resulting in 40 snapshots per simulation. Finally, we note that the initial conditions at lower resolutions ($128^{3}$ and $256^{3}$) have been calculated from the initial $512^{3}$ snapshot of each realization by coarse-graining, i.e. volume-averaging over $2^{3}$ and $4^{3}$ cells, respectively. We choose this approach to minimize the differences in the initial conditions between the individual configurations of a particular realization. In addition to this, all simulations, including the LES, first decay for $0.2T$ without model before the actual $2T$ decay that we follow and analyze. This is done in order to obtain converged spectra at a given resolution, because the process of coarse-graining leaves excess energy at the smallest scales and the interaction between model and excess energy is unknown. Moreover, resolution dependent quantities, e.g. magnetic energy or vorticity (see next section), relax to their intrinsic value in this transient-decay.

III Results

III.1 Energy spectra

Figure 1 shows the kinetic and magnetic energy spectra initially at $t=0T$ and after two turnover times at $t=2T$.

Initially, there is basically no variation of the individual spectra between different configurations. However, the difference in resolution is clearly visible. The highest resolution ensemble (ILES-512) exhibits the most extended power-law range in the kinetic energy spectrum. Accordingly the wavenumber $k$ where the spectrum drops due to numerical dissipation is shifted towards larger scales with decreasing resolution (ILES-256 and ILES-128). This also confirms our approach of removing coarse-graining artifacts in the initial snapshots by the initial transient-decay. We verified that the smallest scales are statistically stationary in the following evolution. The small peak still visible around $k=2$ is due to the original characteristic driving scale of the initial forced simulation.

After two turnover times, the differences between resolutions remain the most striking feature in the spectra with respect to the wavenumbers where the spectrum wears off. At the lowest resolution, the differences between the ILES-128 and the different SGS models are subtle. There exists virtually no difference between the ILES-128 and the grid-filtered eddy-viscosity (EV-GS) and nonlinear (NL-GS) model or the scale-similarity runs (SS-F3 and SS-F5) - both in the kinetic and in the magnetic spectrum. The SGS runs of the eddy-viscosity model with explicit filtering lead to a slightly stronger dissipative behavior. This is expressed by a marginal reduction of energy on the smallest scales $k\gtrsim 30$ and more pronounced ($\approx 40\%$) for the larger filter width (EV-F5) than for the smaller filter width ($\approx 20\%$ - EV-F3). The nonlinear model (NL-F3 and NL-F5) seems to have a very similar dissipative behavior to the eddy-viscosity model given that the spectra coincide for the same filter.

Finally, the simulations are still approximately isotropic after the two turnover times as measured via the generalized Shebalin angle $\theta$ Shebalin et al. (2009); Oughton et al. (1994) defined by

[TABLE]

with $k_{\perp}^{2}=k_{y}^{2}+k_{z}^{2}$. A fully isotropic spectrum yields $\theta=\tan^{-1}\sqrt{2}\approx 54.7^{\circ}$, whereas a fully anisotropic spectrum with all energy in the perpendicular modes yields $\theta=90^{\circ}$. Here, for all simulations, i.e. independent of resolution, SGS model and filter, the generalized Shebalin angle changes from $\theta\approx(55\pm 1)^{\circ}$ at $t=0T$ to $\theta\approx(57\pm 1)^{\circ}$ at $t=2T$ for both the kinetic and magnetic energy spectrum. Thus, the weak mean field is not expected to have a significant influence on the flow over the free decay.

III.2 Temporal evolution of mean quantities

The evolution over time of the spatially averaged kinetic energy, magnetic energy, vorticity magnitude, current density magnitude and dilatation magnitude is shown in Fig. 2. Overall, all quantities smoothly decay over the two turnover times as expected. The panels show a similar behavior of the SGS models as observed in the energy spectra. However, there are subtle differences.

The evolution of the kinetic energy, $\rho\bm{u}^{2}/2$, is converged with respect to resolution and SGS model. In contrast to this, the magnetic energy, $\bm{B}^{2}/2$, shows a clear separation with resolution. The decreased turbulence intensity or effective Reynolds number at lower resolutions cannot sustain the original magnetic field strength of the driven simulation conducted at higher resolution. Thus, the differences in the initial values at $t=0T$ can be attributed to the removal of coarse-graining artifacts (here, the excess magnetic energy for a given resolution) by the transient-decay. When removing the resolution effects, e.g. by normalizing each ensemble to its initial value, all configurations but one (EV-F5) collapse to a converged evolution. The eddy-viscosity model with the largest explicit filter shows a $\approx 10\%$ decrease in magnetic energy indicating increased dissipation. However, in contrast to the energy spectra, here not only the small scales are affected by the model, but a back-reaction onto the largest scales has taken place.

The derived quantities, i.e. the vorticity magnitude $\left|\nabla\times\bm{u}\right|$, the current density magnitude $\left|\nabla\times\bm{B}\right|$ and the dilatation magnitude $\left|\nabla\cdot\bm{u}\right|$, exhibit comparable behavior. In the raw data, resolution effects dominate and resolving smaller spatial scales leads to larger values. Contrary to the evolution of the magnetic energy, this effect does not vanish when all configurations are normalized and a lower resolution results in a slightly increased decay rate. The same four SGS models (EV-F3, EV-F5, NL-F3 and NL-F5) as in the energy spectra now separate from the bulk (EV-GS, NL-GS, SS-F3 and SS-F5), which is indistinguishable from the ILES-128 ensemble. The vorticity and current magnitudes are reduced by $5\%$ (NL-F3), $6\%$ (EV-F3), $10\%$ (NL-F5) and $12\%$ (EV-F5) indicating very similar behavior with respect to filter width between eddy-viscosity and nonlinear model. In contrast, the dilatation magnitude is only reduced by $5\%$ (F3) and $6\%$ (F5) for the NL models, but by $8\%$ (F3) and $16\%$ (F5) for the EV models. The latter indicates that the eddy-viscosity model is more isotropic than the nonlinear one in agreement with their functional form.

III.3 Higher-order statistics

After having described the evolution of mean quantities, we now consider the temporal evolution of the higher order moments of the distributions. They provide insight into the tails of the distributions, which are crucial in the characterization and understanding of, for instance, the intermittency of turbulence. In general, the variances of all quantities (kinetic and magnetic energy, and vorticity, current and dilatation magnitudes) posses the same characteristics as their means in the previous section, i.e. an overall decay is observed with similar ensemble variations and configuration separations.

The next higher order moments we consider are the skewness

[TABLE]

with standard deviation $\sigma$ and the (Fisher) kurtosis

[TABLE]

The kinetic and magnetic energy skewness and kurtosis do not discriminate between the different models as the ensemble variations for each configuration are larger than the differences between the configurations. This picture changes when looking at the higher order moments of the derived quantities. The temporal evolution of the skewness and kurtosis of the vorticity, dilatation and current magnitude are very similar (with respect to the qualitative evolution of the medians and interquartile ranges) as shown in Fig. 3. Thus, we focus on the magnitude of the current density for a quantitative discussion.

Firstly, both skewness and kurtosis are resolution dependent. A lower resolution increases the skewness by $\approx 5\%$ (ILES-256 vs ILES-512) and $\approx 10\%$ (ILES-128 vs ILES-512), and the kurtosis by $\approx 10\%$ and $\approx 25\%$, respectively.

All eddy-viscosity and scale-similarity models follow this trend. They evolve virtually identically to the ILES-128 ensemble. Also the nonlinear model based on grid-scale quantities (NL-GS) does not have a measurable impact on the results. However, the explicitly filtered nonlinear models (NL-F3 and NL-F5) clearly improve over the ensemble without model. Their evolution is consistent with the higher resolution (ILES-256) results. Moreover, the differences between using a three-point stencil and a five-point stencil are negligible indicating a converged result.

To further illustrate this we show the instantaneous probability density function (PDF) of the current magnitude at $t=1T$ in Fig. 4.

The top panel illustrates the raw PDFs. The ensembles at different resolutions are clearly identified by an overall shift on the x-axis. This corresponds to the decrease of the mean with resolution, as described in section III.2. The differences in the higher order statistics of the different configurations are already hinted at in the insets. A pure shift would equally affect the left and right hand side tails. This is observed in the PDFs of EV-F3 and EV-F5, which are both shifted to the left in comparison to the implicit LES configuration. In contrast to this, the PDFs of NL-F3 and NL-F5 only exhibit a shift in the right tail and coincide with the ILES-128 in the left tail indicating a change in the shape of the PDF.

This difference is evident in the PDFs of the normalized current in the bottom panel of Fig. 4, where the individual PDFs have been normalized by the respective mean values. The three features previously identified in the temporal evolution of the skewness and kurtosis, i.e. the resolution differences (ILES-128 vs ILES-256 vs ILES-512), the insensitivity of the EV and SS models, and the improvement by the nonlinear model, are clearly present. In fact, the ensemble distributions of the nonlinear model (NL-F3 and NL-F5) at a resolution of $128^{3}$ match the distribution of the implicit LES at a resolution of $256^{3}$ demonstrating a clear enhancement. Finally, we emphasize again that the results of the current density presented in this subsection are qualitative identical to the ones obtained for the vorticity and dilatation magnitude, i.e. the explicitly filtered nonlinear models match the higher-resolution ILES.

III.4 Structure functions

In order to gain further insight into the flow we analyze the structure functions (e.g. Biskamp, 2003). In particular, we look at the longitudinal velocity structure functions of order $p$

[TABLE]

which are given by the moments of the velocity increments along the direction of separation $\bm{l}$ assuming homogeneity and isotropy. Structure functions are related to the correlation functions and the energy spectrum. Moreover, they exhibit scaling behavior in the inertial range $S^{p}(l)\propto l^{\zeta_{p}}$ so that scaling exponents $\zeta_{p}$ can be determined.

Figure 5 illustrates the second order longitudinal structure function $S^{2}_{\parallel}$ for all configurations of one arbitrary realization after one turnover time of free decay. All structure functions have been calculated based on $10^{10}$ randomly chosen pair of points. The convergence has been verified by comparing the results with the ones obtained by using twice the amount of points for one particular snapshot. Two important features can be observed. First, the structure functions of all configurations, i.e. independent of resolution and presence of an SGS model, collapse (on top of each other) on scales $\gtrsim 30\Delta_{x}$. On smaller scales, the differences with respect to resolution are more pronounced. This is expected since the increasing numerical dissipation with decreasing resolution leads to a decrease of variations in the velocity field on the small scales. Again, the grid-filtered EV and NL LES, and the scale-similarity runs are indistinguishable from the ILES run. The increased dissipation of the EV and NL model already observed in the spectra and mean quantities is also visible here in the slightly reduced variations on the smallest scales. Second, no clear power-law range can be identified in any of the configurations which can be attributed to the limited resolution, which for these simulations indicates a too low Reynolds number.

For this reason, we make use of the concept of extended self-similarity (ESS) stating that the scaling behavior with corresponding scaling exponents can be recovered by relating structure functions to each other. While originally discovered in hydrodynamics (Benzi et al., 1993), this concept works remarkably well in MHD, too. For example, in Fig. 6 we plot $S^{2}_{\parallel}$ versus $S^{3}_{\parallel}$ (of the same snapshot as in Fig. 5). A power-law behavior for all configurations is clearly visible. Moreover, the scaling exponents in this representations are by construction identical to the original ones. Thus, we determine the individual exponents in this representation by nonlinear least-square minimization using the lmfit package (Newville et al., 2014). With these exponents, we now continue our analysis in two directions: reevaluation of the structure functions versus separation distance, and scaling behavior with increasing order $p$.

Figure 7 shows the median (over all 19 realization) second order structure function of all configurations compensated by the corresponding third-order structure function and scaling exponent, i.e. $S^{2}_{\parallel}/(S^{3}_{\parallel})^{\zeta_{2}}$, versus distance $l$. The plot illustrates where and to what extent the power-law scaling is found in the non-normalized data. Approximate power-law scaling is observed for all configurations on scales $\gtrsim 20\Delta_{x}$. Below $20\Delta_{x}$ the individual configurations start to deviate from ideal scaling. The deviations grow towards smaller scales for all configurations, however, to different degrees for different configurations. As expected, the highest resolution runs (ILES-512) exhibit the least deviation (at most $6\%$ on the smallest scale), followed by the intermediate resolution runs (ILES-256) with $\approx 7\%$ on the smallest scale. At the lowest resolution, the no-model (ILES-128), grid-scale filtered SGS (EV-GS and NL-GS) and scale-similarity (SS-F3 and SS-F5) runs show the strongest deviation, $\approx 15\%$. The explicitly filtered eddy-viscosity (EV-F3 and EV-F5) and nonlinear models (NL-F3 and NL-F5) display an improved behavior over the other low resolution runs. While the two EV models deviate by $10\%$, the nonlinear models deviate $8$-$9\%$, reaching almost the performance of the intermediate resolution runs. Qualitatively, the same behavior observed for the structure functions of order $p=2$ is also observed for the structure functions of order $p=1$ and higher orders.

Finally, we analyze how the scaling exponents evolve with order $p$ depending on resolution and presence of an SGS model. Figure 8 illustrates the median coefficients $\zeta_{o}$ over all 19 realizations up to order $p=6$. The coefficients have been determined based on nonlinear least-square fitting and employing extended self-similarity. For higher orders the ESS does not provide robust exponents any more. The left panel of Fig. 8 shows the influence of decreasing resolution on the exponents. While the exponents up to $p=4$ are virtually identical, a lower resolution (ILES-128 and ILES-256) leads to a slight overestimation of $\zeta_{5}$ ($2\%$) and $\zeta_{6}$ ($3\%$) in comparison to ILES-512. For reference, we also plot the exponents as derived by She and Leveque (She and Leveque, 1994) under the assumption that the most singular dissipative structures are filaments

[TABLE]

The reference run ILES-512 fits the prediction remarkably well with a deviation of only $8\permil$ at the highest order $p=6$. In general, the different SGS models and filtering procedures do not have a measurable influence on the scaling behavior. All behave like the ILES-128 yielding slightly overestimated exponents at high order, which can be attributed entirely to the low resolution.

III.5 Computational efficiency

Finally, we compare the additional computational costs incurred with the calculation of the different SGS models and the filtering. Given that all LES were conducted at a resolution of $128^{3}$ grid points, we compare the relative overhead over the no-model ensemble at the same resolution (ILES-128). Table 3 lists the mean ratios of the time per individual cycle at the fluid level, i.e. other factors such as inter-process communication are not included.

The time per cycle increases for all SGS models when compared to the ILES-128 as expected. Furthermore, two general trends are visible.

First, the computational costs increase with increasing filter width. For example, the nonlinear model with grid-scale quantities increases the time per cycle by a factor of $\approx 1.18$ (NL-GS). Explicit filtering introduces additional computations and is thus even more expensive, i.e. a factor of $\approx 1.47$ for NL-F3 and of $\approx 2.41$ for NL-F5, respectively. The unproportional increase in computational costs between F3 and F5 is easily explained by the unfavorable memory access in the filtering procedure. For F5 the filter is build upon a stencil involving $5^{3}=125$ points resulting in many accesses to non-contiguous memory and thus cache-misses.

Second, the eddy-viscosity and nonlinear model introduce a similar overhead, with EV being a few percent cheaper than NL, while the scale-similarity model is about twice as expensive as the other two models. The latter is attributed to the additional (comparatively expensive, explicit) filter operations. NL and EV only require filter operations on the 7 primary quantities ($\overline{\rho},\widetilde{\bm{u}}$ and $\overline{\bm{B}}$). The scale-similarity models also needs all filtered mixed quantities (\mathop{\vbox{\halign{#\cr\kern 3.0pt\cr\bracelu\leaders\vrule height=0.39998pt,depth=0.0pt\hfill\braceru$\crcr\kern 3.0pt\nointerlineskip\cr$\hfil\displaystyle{\widetilde{u_{i}}\widetilde{u_{j}}}\hfil\crcr}}}\limits, \mathop{\vbox{\halign{#\cr\kern 3.0pt\cr\bracelu\leaders\vrule height=0.39998pt,depth=0.0pt\hfill\braceru$\crcr\kern 3.0pt\nointerlineskip\cr$\hfil\displaystyle{\overline{B_{i}}\overline{B_{j}}}\hfil\crcr}}}\limits, and \mathop{\vbox{\halign{#\cr\kern 3.0pt\cr\bracelu\leaders\vrule height=0.39998pt,depth=0.0pt\hfill\braceru$\crcr\kern 3.0pt\nointerlineskip\cr$\hfil\displaystyle{\widetilde{\bm{u}}\times\overline{\bm{B}}}\hfil\crcr}}}\limits) which involves 15 additional filter operations in total.

Finally, we also tested how the computational costs increase for a no-model run at the next higher resolution (ILES-256). At the level of a single fluid cycle the time increases by factor of $\approx 9$. However, this does not yet take into account that the timestep is also reduced by a factor of $\gtrsim 2$ at $256^{3}$ versus $128^{3}$. Thus, the total time required to reach a certain state in the simulation is effectively increased by a factor of $\gtrsim 18$.

IV Discussion

One of the most striking results from the analysis in the last section is that models calculated from quantities at the grid-scale, i.e. EV-GS and NL-GS, and scale-similarity models (SS-F3 and SS-F5) have no measurable impact on the statistics of the flow. The results for the grid-scale based models are in agreement with findings for finite-difference schemes (Ghosal, 1996; Chow and Moin, 2003) and for shock-capturing methods (Garnier et al., 1999), i.e. numerics dominate over (eddy-viscosity type) SGS models when no explicit filtering is applied. However this does not explain the results for the scale-similarity closures, which employ an explicit filter. A possible explanation for the absence of any observable effect (apart from a very short transient behavior, e.g. in the kurtosis and skewness of the current at $t<0.1T$ as visible in Fig. 3) is that the filter separation is still too small. Physically, this translates to the statement that the modeling assumption of self-similar turbulence is not fulfilled on these numerically strongly damped scales.

Another observation concerns the convergence with filter width for explicitly filtered eddy-viscosity (EV-F3 and EV-F5) and nonlinear (NL-F3 NL-F5) models. The energy spectra and mean quantities exhibit a small dependency on the filter width indicating increased dissipative behavior with larger $\Delta$. However, this is secondary from a practical point of view because a smaller explicit filter is desirable for two reasons. First, higher order statistics, e.g. skewness, kurtosis and normalized structure functions, show approximately converged results for F3 and F5. Thus, the improvements over the ILES-128 can already be achieved with the smaller explicit filter and the nonlinear model (NL-F3) while possessing a smaller intrinsic dissipation and being computationally more efficient than the F5 counterpart. Second, larger explicit filters in their current form are impractical for actual LES in any case.

Due to the unfavorable memory access in the filtering, the computational cost grows exponentially with increasing filter width. In addition, more and more ghost zones are required increasing the costs even further. Building multi-dimensional filters based on the simultaneous application of one dimensional ones rather than the sequential application could be a potential way out. Even though they are found to be slightly less accurate (Sagaut and Grohens, 1999), their multi-dimensional stencil size is decreased dramatically, i.e. from $N^{3}$ to $3(N-1)+1$ supporting points in three dimensions with $N$ being the number of points for the 1-d filter. Alternatively, the filtering could be realized in spectral space. While the process of filtering itself is then reduced to a simple local multiplication, additional complexity independent of the filter width, is introduced by the transformations between real and spectral space.

Our current filtering framework could also be further optimized to reduce the computational overhead of the filtering. For example, cache misses would be partly avoided by using fixed, compiled-in stencils rather than dynamic ones in each cycle. Independently of this, in practice the estimated SGS modeling overhead in a simulation is rather conservative. The total wallclock time always depends on additional factors than the time spent at the pure fluid level, most notably inter-process communication. Moreover, additional physics present in the simulation such as gravity or chemistry can reduce the relative overhead introduced by an SGS model even further.

All models and filters lead to stable simulations and we did not employ any explicitly regularization. While this comes as no surprise for the eddy-viscosity type models, which are only capable of transferring energy down-scale, other groups, e.g. (Winckelmans et al., 2001; Müller and Carati, 2002), typically find that regularization is required for scale-similarity and nonlinear type models. These models also allow for up-scale energy transfer and are thus capable of seeding numerical instabilities when this inverse transfer is not controlled. Most frequently, both type of models are therefore supplemented with an additional eddy-viscosity type term which successfully stabilizes the simulations. However, this only concerns (non shock-capturing) finite-volume, finite difference or (pseudo-) spectral schemes. In our case of a shock-capturing finite-volume scheme, the inherent numerical dissipation acts as an effective eddy-viscosity model (thus the term implicit LES, see e.g. (Grinstein et al., 2007)) and evidently provides sufficient regularization for stable simulations.

It should be noted the present study only analyzes turbulence with negligible cross-helicity. In order to account for cross-helicity effects such as changing cascade dynamics both structural models are expected to be applicable as presented. While the nonlinear model makes no assumptions on the nature of the flow Vlaykov et al. (2016), the scale-similarity model implicitly accounts for changing dynamics Chernyshov et al. (2012). The functional eddy-viscosity model misses this feature. For this reason, extensions have been proposed to explicitly capture unresolved cross-helicity effects Yokoi (2013). Similarly, all models presented concern fully collisional plasmas. Kinetic effects from low collisionality plasmas such as anisotropic thermal conduction even in the presence of weak mean fields or heating and nonthermal particle acceleration from magnetic reconnection are not captured by the presented models Miesch et al. (2015).

Finally, the resolution of the LES ($128^{3}$) in this work prohibits the study of more detailed physical effects in MHD turbulence itself Zhou and Oughton (2011) and the corresponding properties of the SGS models. For example, in order to analyze the locality of interactions Zhou (2010) a much higher resolution is required so that the effects of numerical dissipation can be clearly distinguished from physical ones. In addition, higher resolution LES would also allow to study the properties of the SGS models with respect to how they affect the straining (i.e. the distortion of small scale vortices by large scale motions) and sweeping (i.e. advection of small scale vortices by the large scale flow with negligible distortions) in MHD turbulence Zhou et al. (2004). However, we expect that with increasing resolution the importance of a pure SGS turbulence model decreases as more and more features are naturally resolved.

V Conclusions and outlook

In this paper we analyzed the free decay of homogeneous, isotropic, supersonic MHD turbulence with different SGS models and without explicit model on various grid resolutions. SGS models are typically introduced to LES in order to incorporate unaccounted effects from below the grid-scale, and to improve the quality of the simulation at lower computational cost. We measured the quality of the SGS models by their capability to reproduce the results of a reference simulation at higher resolution. The reference quantities included the energy spectra, the evolution of different statistical moments of the kinetic and magnetic energies, the vorticity, the current density and the dilatation magnitudes. In total, we compared three SGS models: eddy-viscosity, scale-similarity and nonlinear. Additionally, we evaluated the influence of using implicit-, grid-filtered quantities versus explicitly filtered quantities to calculate the model terms. We analyzed an ensemble of 19 different initial conditions for each configuration as temporary, transient fluctuations can easily dominate individual simulations.

We find that the simulations employing a grid-filtered eddy-viscosity (EV-GS) or nonlinear (NL-GS) model, or a scale-similarity model with the tested explicit filters (SS-F3 and SS-F5) produce results that are indistinguishable from an implicit LES, i.e. without an explicit model, at the same resolution. Moreover, we find that the eddy-viscosity and nonlinear models with the two tested explicit filter widths, i.e. with filter widths of $2.71\Delta_{x}$ (EV-F3 and NL-F3) and $4.75\Delta_{x}$ (EV-F5 and NL-F5), introduce little additional dissipation on the smallest represented scales e.g. in the energy spectra or the evolution of the mean quantities. Finally, the nonlinear model (NL-F3 and NL-F5) improves higher order statistics of small-scale dependent quantities, such as the kurtosis and skewness of the current density, dilatation and vorticity. For these quantities, the results of an ILES at doubled resolution (in each dimension) can be achieved while introducing only a small computational overhead — less than factor of $1.5$ (versus $\approx 16$ for the higher-resolution ILES). This similarly applies to the normalized structure functions and is independent of the explicit filter width.

Based on these results we conclude that an explicit filter is required in order to obtain a measurable impact of an SGS model for shock-capturing finite-volume schemes of second order. In how far this conclusion holds for schemes of higher order and more dynamic versions of the SGS models, e.g. with dynamic coefficients, is yet to be seen and subject to future work. Furthermore, additional dissipation for the explicitly filtered models is not required as numerical dissipation proves to be sufficient (if not too high). Thus, the introduction of an eddy-viscosity model in these schemes is unnecessary. However, as the nonlinear model improves higher order statistics, it would be desirable to remove the unnecessary dissipation to improve the lower-order statistics as well. This kind of regularization is also subject of future work. Finally, the nonlinear model in its current version can readily be used with a small explicit filter in situation where higher order statistics are important at little extra cost. The associated code will be made publicly available together with the publication.

Acknowledgements.

PG acknowledges financial support by the International Max Planck Research School for Solar System Science at the University of Göttingen. DV acknowledge research funding by the Max-Planck-Institut für Dynamik und Selbstorganisation. DRGS thanks for funding through Fondecyt regular (project code 1161247), through the ”Concurso Proyectos Internacionales de Investigación, Convocatoria 2015” (project code PII20150171), and from the Chilean BASAL Centro de Excelencia en Astrofísica y Tecnologás Afines (CATA) grant PFB-06/2007. The Enzo simulations were performed and analyzed with the HLRN-III facilities of the North-German Supercomputing Alliance under grant nip00037.