Fitting of extended subgrid scale models in compressible turbulent MHD
Abstract
The use of Large Eddy Simulations in compressible, turbulent MHD is often limited to the implementation of purely dissipative subgrid scale models. While such models work well in the hydrodynamic case due to the universality of the spectrum, they do not fully describe the complex dynamics of MHD, where the transfer of energy between internal, kinetic and magnetic energies at small scales are less trivial. For this reason, a subgrid scale model based on the gradients of the fields entering in each nonlinear term of the equations has already been proposed and studied in the literature, for the momentum and induction equations. Here, we consider the full set of compressible, ideal MHD equations, with an ideal gas equation of state, and we proposed a generalization of the gradient model, including the energy equation. We focus on the residuals coming from the whole set of equations, by filtering accurate highresolution simulations of the turbulent Newtonian KelvinHelmholtz instability in a periodic box. We employ the same highresolution shock capturing methods typically used in relativistic MHD, applicable in particular to neutron star mergers. The apriori test, i.e. the fit between the subfilter residuals and the model, allows us to confirm that the gradient model outperforms any other, in terms of accuracy of the fit and small deviations of the bestfit precoefficient from the expected value. Such results are validated for 2D and 3D, for a range of different problems, and are shown, for the first time, to hold also for the energy evolution equation. This paper is the first step, based on a solid theoretical and numerical basis, towards the nearfuture extension of the subgrid scale gradient model to the relativistic MHD, and a future implementation in a full General Relativity LES.
I Introduction
In the context of computational fluid dynamics, direct numerical simulations (DNS) are employed to resolve the equations in all the dynamically relevant scales, including the dissipative ones. This means that the numerical grid size has to be smaller than the dissipative scale. Since the computational cost of DNSs scales as (Sagaut, Deck, and Terracol, 2006) (where is the Reynolds number), for many problems (including convection zones of stars, planetary atmospheres, accretion disks, and industrial applications such as gas turbines, steady turbulent flows…) it is unfeasible to have a small enough grid size and a large enough domain at the same time, thus representing all the scales of interest. Even though in the near future the increasing computation power will improve the achievable resolution, thus making DNSs feasible for more scenarios, currently there are many cases where we are still very far from it. In astrophysics, one of these cases, which is of particular interest due to the recent discoveries of gravitational waves, is the binary neutron star merger.
The two most widely used numerical approaches to overcome the lack of resolution in turbulent flows are the Reynolds Average Numerical Simulations and the Large Eddy Simulations (LESs). For both cases, the idea is to directly evolve the quantities on the resolved scales, i.e. from the integral scales (the largest ones, corresponding to dominion size or to the typical scale where energy is injected), down to the grid scales, and to describe the effect of the smaller ones through a suitable subgrid scale (SGS) model.
While the purely hydrodynamic turbulence dynamics is relatively well understood due to the universality of the spectral properties, the dynamics becomes more problemdependent when magnetic fields are included. Magnetic fields tend to break the usually assumed isotropy and selfsimilarity, and prevent one from defining universal laws, valid for different scenarios. Magnetohydrodynamic (MHD) turbulence includes a nonlinear dynamics, with transfer between magnetic and kinetic energy (dynamo mechanism, when the magnetic energy rises), and, at the same time, between different scales. Such transfers vary locally and in time, so that a statistical evaluation of them represents essential information. In general, the smallest scales are the most interesting ones in terms of interchange of kinetic and magnetic energy.
As a consequence, it is crucial to properly model the unresolved scale dynamics, by means of a SGS modeling, that is incorporated as additional extra terms in the MHD equations. Due to the variety of turbulent scenarios, LESs with SGS models represent an active field of interdisciplinary research, including engineering, plasma physics and astrophysics (i.e., see for instance a review about the use of LESs in astrophysicsMiesch et al. (2015), which is also our main interest). A popular SGS model is the dissipative one, introduced by Smagorinsky half a century ago in the context of hydrodynamic direct cascadesSmagorinsky (1963) to mimic the fluid’s viscosity. Other SGS models have also been proposed relying mainly on dimensional analysis or physical assumptions like the dynamo mechanism.
Another branch of SGS models (including the multiscale and multivariational ones)Hughes et al. (2001) relies on selfsimilarity arguments, where basically the loss of information is estimated considering the average spatial correlations of the evolved fields at different scales, and extrapolating such information to infer the fields in the SGSs. These models are used especially in finiteelements schemes.
On the other hand, the gradient modelMüller and Carati (2002) is based on a mathematical formulation which relates the finite resolution to an effective filtering of the equations. By expanding in Taylor series the nonlinear terms in any equation, one can approximate the SGS residuals with a functional form which only relies on mathematical assumptions. The model does not contain any arbitrary physical assumption, and just extrapolates to the SGS the trends seen at the smallest resolved scales.
In general, two steps can be taken in order to assess and compare the performance of any SGS model:

Apriori fitting. This test require the postprocess filtering of the fields obtained with a relatively highresolution simulation at a given time. The residuals coming from the filtering (i.e., the nonlinear information contained between the filter size and the original resolution ) can be fit to a given SGS model, where the only free parameter is a precoefficient. The procedure can be easily repeated for different times, resolutions, filter sizes, initial conditions, in order to test the applicability of the tested SGS model. This procedure allows to identify the bestfitting functional form of the SGS model, but does not allow to evaluate the effect of the inclusion of the SGS in a simulation.

LES comparison, including different SGS models. Ideally, a lowresolution LES with a good SGS model should resemble a LES with a higher resolution. The LES comparison is a more demanding test, since it involves the dynamical feedback given by the implemented SGS terms. This test is much more informative and is crucial to evaluate the physical meaning of including a certain SGS tensor in a LES. However, the exploration of parameters is much more computationally expensive.
In this paper, we will perform the first step: an extensive apriori assessment of the most popular explicit local SGS models present in the literature, leaving a battery of LES comparison tests for a forthcoming paper, which will obviously rely on the results found with the first approach. Works in this sense have been performed in the context of either incompressibleKessar, Balarac, and Plunian (2016) or compressible Grete et al. (2016, 2017a, 2017b); Grete (2017) MHD. These studies are among the most detailed in comparing different SGS models, employing either forced or decaying turbulence. Although they consider the momentum and induction MHD equations, the internal energy evolution is not included, since they impose either incompressibility or a polytropic/isothermal equation of state (EoS) as a closure Grete et al. (2017b). Here we extend these works by considering instead the compressible MHD equations with an ideal EoS.
We run several bounding box highresolution simulations of the KelvinHelmholtz instability (KHI), which triggers the gradual development of turbulence. We do not include any forcing: we focus on dynamical, decaying turbulence, since we aim at applying them to simulations of astrophysical transient phenomena. We perform apriori tests, where the filtered residuals are compared to what we would obtain with different SGS models applied to the information available at the filter resolution. We validate both the method and the SGS models, identifying the bestfit coefficients for a variety of initial data, resolutions and filters.
In § II we review and extend the formal description of a LES, seen as a consequence of separating scales by applying a filter to the solution. We describe in some detail the filtered full MHD equations in § III. In § IV we introduce different SGS models, with special emphasis on the gradient model. In § V we describe the platform and the numerical methods used in this work (in Appendix A we describe the validation of the code with a set of benchmark tests in 2D). In § VI we present the main highlights of the results coming from 3D KHI simulations, showing a variety of cases. Conclusions are drawn in § VII.
Ii Equivalent filtering in Large Eddy Simulations
In a LES, the smallest simulated scales (i.e., the ones corresponding to the grid size ) are not sufficient to capture entirely all the relevant dynamics, which might occur in the SGS scales. Since the SGS part of the solution is unknown by definition, it can only be either neglected, leading to the socalled implicit LES,gri (2007) where a lowpass filterSun and Domaradzki (2018) or the internal numerical dissipation (used for instance in solar physicsGhizaru, Charbonneau, and Smolarkiewicz (2010)) act as an intrinsic SGS model, or explicitly modeled through the information which is available, i.e. the one of the resolved scales. This explicit modeling, on which we focus, can be based either on mathematical or physical arguments.
Mathematically, the effect of finite resolution can be thought as equivalent to a lowpass spatial filter, where the filter size, , is given by the grid cell, . This separates, for a continuous field , the smoothed (resolved) part, from the subfilter scale (SFS) fluctuations (or residuals), respectively:
(1) 
Indicating the filter kernel with , the filtering operator over a field can be written as
(2) 
An homogeneous isotropic lowpass filter is independent on the direction (i.e., ) and only smooths out fluctuations on length scales smaller than the filter size , leaving unchanged the variations of the solution at larger length scales. In addition, the filter operator is linear and commutes with spatial derivatives. Generically, it can be written for any dimension as
(3) 
where is just the onedimensional kernel function.
The simplest lowpass filter is the mean value in a cubic domain with size in each Cartesian direction , described by the normalized kernel
(4) 
Despite the appealing simplicity of the box filter, which makes it very useful to perform numerical calculations, we will see below that it is not suitable for analytical calculations involving its derivatives, since they are not continuous. Therefore, at a formal level, it is more practical to introduce the normalized Gaussian kernel, which in the space domain can be written as
(5) 
where defines the effective filtering width. Besides having the same zeroth and first moments, Gaussian and box filters have the same second moment if we set :
(6) 
The filter is a useful mathematical tool to analyze not only the different scales on the solution, but also the structure of the evolution equations, allowing us to distinguish between the resolved or filtered fields and the unresolved ones on the SFS.
In order to see the effect of this scaleseparation in a simple application, let us consider the wellknown Burgers equation, a longstanding field of research in LES.D. Love (1980) In the 1D case, given the velocity field , the equation reads:
(7) 
The equation for the filtered fields can be obtained by applying the filtering operator to the equations. Since the filtering operator commutes with both spatial and time derivatives, we have
(8) 
Notice that the filtered equation is equivalent to the original one except by a term proportional to which represents the SFS residuals, that is, the loss of smallscale information due to the filtering process related to the non linear terms. If in the region over which we average has a definite sign (which is usually the case except close to the shock), then and is negativedefinite, so that represents the kinetic energy hidden in the SFS.
Regardless of the specific equations considered, these new terms always appear due to the noncommutativity of the filtering operator with the nonlinear terms of the equations. In a LES, which effectively only evolves and cannot simulate the unresolved scales, the explicit expression of is apriori unknown, and needs to be written as a function of the filtered fields in order to close the system of equations.
There are several ways of doing that, based mostly either on physical arguments or on an expected selfsimilarity of the solution. For instance, a popular, historical approach to fluid dynamics SGSsSmagorinsky (1963) is the setup of an artificial viscosity , namely
(9) 
The viscosity parameter is often taken to be proportional to due to dimensional reasons and in order to ensure that this term vanishes in the continuum limit, thus guaranteeing numerical convergence to the continuum (or DNS) solution. The numerical value of the proportionality coefficient can be fixed by hand or estimated by means of dynamical procedures which assume selfsimilarity (thus similar to a multiscale/multivariational approach).Germano et al. (1991) Regardless of the value of , this functional form of is equivalent to the viscous Burger’s equation, where the viscosity is the numerical SGS model.
However, in many cases the physics involved does not consist only in dissipation, but can involve, for instance, inverse cascade, and scaledependent transfers of energy between the fluid kinetic energy and the magnetic one. Therefore, in those cases any result coming from the LES will be biased by the apriori choice of a given SGS physics.
If one wants to avoid introducing arbitrary physical assumptions, it is possible to compute these terms relying only on mathematical arguments by considering the analytical Taylor expansion of the SFS termsLeonard (1975); Müller and Carati (2002). The main hypothesis is that the discretization of the equations over a finitesize grid is approximated by a filtering operator acting on the equations, with a Gaussian kernel equivalent to a box filter with width , eq. (5). Its dimensional Fourier transformed function is
(10) 
The main idea is to compute an approximation of the inverse filtering operator based on gradient expansion of the filter kernel , that is, an approximation of the inverse Fourier transform of . The first step is to perform a Taylor expansion of the transformed function and its inverse in terms of the filter scale, that is,
(11)  
(12) 
Considering the expansions to the fields and , the application of the inverse Fourier transformation yields to an infinite series representation of the filter operator and its inverse in terms of gradient operators acting on the fields, namelyVlaykov et al. (2016)
(13)  
(14) 
These expressions are absolutely convergent and formally accurate at all orders, since the Gaussian kernel is infinitely differentiable and with unbound support. In fact, it was found that these series converges for all canonical filters.
The unknown components of the SFS tensors usually have the form and . Applying the relations in Eq. (14) results in a series in terms of the individual fields, namely
(15)  
which is accurate up to . The first order terms scale with the second moment of the filter, eq. (6). With this expansion is easy to model the SFS terms appearing in Eq. (8), namely
(16) 
For a Cartesian grid with step , the prefactor can be set to the equivalent value , according to eq. (6). In principle, corrections to these expression are expected due to the form factor (due to the kernel shape) and the contribution from higher orders. As a matter of fact, if the range of modeled SFSs, , becomes large, then the higherorder terms can be important. In other words, this approximated expression can help as long as we can capture most of the dynamics with our LES, leaving to the SGS the task of mimicking the smallscale contributions.
Regardless on these quite obvious caveats, notice that this prescription is very different from the viscous one given by eq. (9): one arises from physical considerations, while the other just from mathematical ones. After this simple example, we pass to the application of this formalism to the MHD equations.
Iii Filtering the Newtonian compressible MHD equations
iii.1 Compressible MHD equations
We consider the conservative formulation of the ideal Newtonian compressible MHD equations, consisting of the continuity, Euler, induction and energy evolution equations, respectively:
(17)  
where the magnetized fluid is described by the total energy density defined as
(18) 
the mass density , the fluid pressure , the internal energy density , and the velocity and magnetic field vectors, being the standard Kronecker delta. For simplicity, we have set the speed of light and the magnetic permeability . The conserved quantities of the equations , in their discretized versions, are numerically evolved. At each step of a simulation, the internal energy is recovered by inverting eq. (18).
The pressure, required to close the system of equations, is defined by an EoS, which in general can read . In this paper, we will present results with the ideal gas EoS, namely
(19) 
where is the ideal gas adiabatic index. However, our formalism is general and can extend to any EoS.
iii.2 Filtered continuity equation and variable inversion
The filtered continuity equation is simply
(20) 
where the Favrefiltered velocities, i.e., weighted by density(Miesch et al., 2015)
(21) 
are the natural choice in compressible MHD, where the momenta are the evolved fields. Therefore, and are simply written as under the filtering operation.
We will use the same notation whenever we indicate an auxiliary field as a function of evolved filtered fields only. Specifically, in our case we can express the internal energy density as a function of the evolved filtered variables, eq. (18), defining
(22) 
Note that there is a difference between what can be numerically evaluated with the resolved fields only, (i.e., combining the filtered evolved fields), and what cannot be, (i.e., filtering the local nonlinear combinations of the variables):
(23) 
which can be interpreted as the kinetic and magnetic energy density hidden in the SFSs. Note that, if such SFS contributions (usually positivedefinite) are neglected in the variable inversion, this will cause a local artificial rise of the evaluated internal energy, .
Similarly, we can define the pressure as a function of all the other evolved variables . Its expression depends on the EoS, as we will explain below.
iii.3 Filtered momentum and induction equations
On the other hand, the filtered momentum and induction equations become, respectively
(24)  
(25) 
where , and indicate the tensors arising from the differences in the non linear products , and , respectively, when the filtering is considered:^{1}^{1}1In literature, different sign conventions are used. We use the same sign convention as Grete’s papers.Grete et al. (2015, 2016); Grete (2017) For instance, in Kessar’s workKessar, Balarac, and Plunian (2016), , are definite with an opposite sign, while has the same sign.
(26)  
(27)  
(28)  
(29) 
Notice that, by construction, and are symmetric tensors, while is antisymmetric and its inclusion is equivalent to having a SFS electric field with components given by . Note that the SFS kinetic and magnetic energies defined in eq. (23) are related to and by
(30)  
(31) 
In the literature, sometimes the functional forms of the traces of and are defined separately from the ones of the offdiagonal components (related to the strain). This is useful as a tight analogy with the standard bulk and strain viscosity, but, since our formulation is general and does not rely on the underlying physics, we prefer to avoid this distinction.
The term associated to pressure is instead purely diagonal and has been only partially considered before. If isothermal EoS or incompressibility are assumed, pressure does not depend on the evolved field, and the only surviving term is the magnetic pressure, given by half the trace of . However, in the general case, the pressure depends nonlinearly on the evolved fields, so that SFS residuals implicitly appear and is generally not zero. In the case of the ideal EoS, such term is proportional to eq. (23), so that:
(32) 
Note that an ideal quasiisothermal gas () can be taken as a shortcut to minimize the spurious increase of pressure.
iii.4 Filtered energy equation
While the momentum and induction equations have been explored in many different works, the energy evolution has usually been neglected (for incompressible and isothermal cases), or only partially considered. In particular, some authors have focused on the SFS kinetic and magnetic energy, studying in detail the interscale transfers between resolved and SFSsKessar, Balarac, and Plunian (2016); Vlaykov et al. (2016); Grete et al. (2017b). Those works aim at analyzing the energy transfer, but not to perform a detailed modeling of the SFS tensors appearing in the evolution equations. Hereafter, we propose two approaches. The first is a conservative one, similar to the previous equations, where only flux terms appear, and is apt for the implementations of SGS models in a LES. In the second one, we considered the already studied interscale energy transfers, appearing as source/sink terms.
Conservative approach (for SGS modeling in LES).
By following a similar approach as for the momentum and induction equations, we can write the energy evolution equation as
(33) 
where the additional vectorial SFS terms on the righthand side arise from the different nonlinear terms in the original equation, as follows:
(34)  
(35)  
(36)  
(37)  
(38) 
We distinguish three fundamental contributions. The first one is given by the SFS pressure advection term (), and the SFS convective term (). Combined, they represent the advection of the SFS enthalpy:
(39) 
where we have exploited the identity to introduce the enthalpy density . In the case of an ideal EoS, this relation is reduced to
(40) 
where the internal energy and pressure as functions of the filtered evolved quantities are given by eq. (22).
The second group of SFS terms in the right hand side of eq. (33) is given by the SFS advection of kinetic () and magnetic energy () densities. Note that, strictly speaking, the latter actually comes for one half from the magnetic energy advection, and for the other half from the magnetic pressure advection, similarly to . In order to keep the elegance of the enthalpy term described above, we prefer to keep the magnetic pressure inside this term.
The last term in the right hand side of eq. (33) is the crosshelicity term, which is partially related to the induction equation tensor, , but cannot be written only as a function of it (or a simple contraction with another filtered quantity).
The main disadvantage of this approach is the presence of triple products of three or four conserved fields in the residuals, which makes the calculation not straightforward. In this paper, we will not focus on this approach, which will be studied in depth in our generalized formulation to be published in a forthcoming paper.
Interscale approach (for analysis and fitting).
Alternatively, we can consider the evolution in time of the different components of the total energy density, . After some algebra,
one can write the final expression
(41)  
where .
This second approach is useful as an analysis tool, in the apriori fitting tests, as in this paper.
iii.5 Final filtered equations
With these considerations, we can summarize the final set of equations to be evolved for the filtered fields are:
(42)  
(43)  
(44)  
(45)  
where
(46)  
(47) 
being
(48)  
(49)  
(50)  
(51)  
(52) 
In the energy equation, we have reorganized the SFS source terms. Due to the (anti)symmetric nature of the SFS tensors, we have that the integration by parts leaves a part of the terms as flux divergences, and another part which is the contraction of only the symmetric part of (i.e., the resolved fluid strain tensor ), and the antisymmetric part of (i.e., the resolved current tensor ). The sink term is composed by the transfers of advected pressure, kinetic and magnetic energy from the resolved scales to the SFS, respectively:
(53)  
(54)  
(55)  
(56) 
In analogy with the viscous and resistive terms in nonideal MHD (proportional to the strain and current tensors, respectively), at each point of the dominion represents the transfer of energy from the resolved scales into the SFSs (positive values) or viceversa (negative).
Note that, with this compact formulation, we need to give the prescriptions for four tensors, eqs. (48)(51) and two scalars, related with pressure, i.e. eq. (54) and the term in eq. (52). The latter is automatically determined by and in the case of an ideal EoS. Once these six SFS tensors (having a total of 20 independent components in 3D) are prescribed, the sources and are as well, since they do not contain additional degrees of freedom.
Hereafter, to avoid confusion, we will indicate with an overline the SFS residuals tensors resulting from either the formal analysis or the postprocess filtering of a numerical solution, and with a simple the analytical models of the SGS tensors that can be implemented in LES, which is the subject of the next section.
Iv SGS models
The explicit expression of the SGS tensors represents the closure of the system of equations in a LES. By definition, in a LES it is impossible to know the exact numerical values or the functional forms of the SFS tensors and sources above, and . The SGS models considered here are local, i.e. depend only on the derivatives of the fields at the smallest resolved scales: we do not consider the multiscale (or selfsimilar) models, where a second filter is performed in order to evaluate the SGS model. We also neglect the structure models, based on different possible combinations of mutual contractions of the four firstderivative tensors (fluid strain, magnetic strain, current, vorticity), maintaining the parity needed by SGS models in the equations.(Speziale, 1991; Kosović, 1997; Lu and Rutland, 2016) In general, they do not provide a good fit (apriori tests) and are also known to be numerically less stable when implemented in a LES. Below we list the ones we explicitly compare in this work.
iv.1 Gradient model (extended formulation)
One can consider the analytical Taylor expansion of the SFS terms appearing in the MHD equations, under the hypothesis of having a filter with a Gaussian kernel. To apply the presented derivation selfconsistently to the compressible SFS stresses, one has to take into account the compressibility effects from the massweighted filtering operator (i.e., the Favre filtering), which is given by
(57) 
Using the previous relations one can write
(58) 
One can solve for , derivate this expression and use the recurrence relation to obtain, at linear order in ,
(59) 
so that the following expression, function of instead of , holds
(60) 
The correction, is zero in the incompressible case. Using the rules described before, the expansions for the SFS terms (48)(50) give, at leading order in :
(61)  
(62)  
(63)  
This formulation, with as in eq. (6), is known as the gradient modelMüller and Carati (2002) (sometimes called generically nonlinear(Grete et al., 2016)) for MHD.
If we consider the dependencies of the total (fluid plus magnetic) pressure for an ideal gas with index , as given by eqs. (23) and (32), then eqs. (51)(52) are expressed by:
(64)  
(65)  
Besides the wellknown magnetic pressure contribution , there is an additional term, related to the fluid pressure , proportional to times the SFS kinetic and magnetic energy. Note again that for an isothermal EoS (), the fluid contribution vanishes and the only term left is the one related to the magnetic pressure, in agreement with previous formalism.Grete et al. (2016) The proper modeling of includes a much more complicated expression. Since in a LES we will use only conservative terms (i.e., no source terms will be modeled), we will avoid to model and study this term in detail.
While the induction and momentum equations obtained had already been explored in depth,Miesch et al. (2015); Grete (2017) the pressure and enthalpyrelated terms in the momentum and energy equations are a novelty, as far as we know.
The gradient model, despite its relative simplicity, strongly differs from the linear models, explained below and largely employed in literature, because there is no apriori dissipative nature. Moreover, since they reproduce the dynamics seen as the smallest resolved scales where the interenergy transfers are more effective, they allow a strong coupling between velocity and magnetic fields.
iv.2 Dissipative model (Smagorinsky)
The basic purely HD model (Smagorinsky model)Smagorinsky (1963) accounts for the viscosity at small scales, introducing an effective viscosity tensor given by
(66) 
where is a factor to be determined, and the proportionality to is taken for dimensional reasons. For the isotropic HD turbulent cases, where the universal laws apply, the value of numerically in spectral methods, being around (Lu and Rutland, 2016; Germano et al., 1991), depending also on the problem. The coefficient can be determined with a dynamical procedure, instead of being fixed, as done for a channel flow.Germano et al. (1991); Oberai and Wanderer (2005)
A subsequent workTheobald, Fox, and Sofia (1994) extended the Newtonian Smagorinsky model to the MHD, introducing the analogous turbulent magnetic diffusivity in the induction equation, proportional to the current density tensor. We will test then the following tensors:
(67)  
(68)  
(69) 
which are the ones commonly used in the literature, properly adapted (factors with introduced on dimensional basis) and extended with the magnetic tensor , that has been introduced only by analogy to the fluid strain, but do not corresponding to any dissipative mechanism. In literature, different normalizations are present, taking trace and offdiagonal components separately, or considering different dependences with . However, they do not show in general a substantial, systematic improvement,Grete et al. (2017a) and we only test the expression above.
iv.3 Crosshelicity model
This model(Müller and Carati, 2002) considers the dissipation of crosshelicity, which is given by two contributions, proportional to the fluid and magnetic strain rate tensors, , and to the product of current and vorticity , so that