Coherent vorticity and current density simulation of three-dimensional magnetohydrodynamic turbulence using orthogonal wavelets

A simulation method to track the time evolution of coherent vorticity and current density, called coherent vorticity and current density simulation (CVCS), is developed for three-dimensional (3D) incompressible magnetohydrodynamic (MHD) turbulence. The vorticity and current density fields are, respectively, decomposed at each time step into two orthogonal components, the coherent and incoherent fields, using an orthogonal wavelet representation. Each of the coherent fields is reconstructed from the wavelet coefficients whose modulus is larger than a threshold, while their incoherent counterparts are obtained from the remaining coefficients. The two threshold values depend on the instantaneous kinetic and magnetic enstrophies. The induced coherent velocity and magnetic fields are computed from the coherent vorticity and current density using the Biot–Savart kernel. In order to compute the flow evolution, one should retain not only the coherent wavelet coefficients but also their neighbors in wavelet space, and the set of those additional coefficients is called the safety zone. CVCS is performed for 3D forced incompressible homogeneous MHD turbulence without mean magnetic field for a magnetic Prandtl number equal to unity and with 2563 grid points. The quality of CVCS is assessed by comparing the results with a direct numerical simulation. It is found that CVCS with the safety zone well preserves the statistical predictability of the turbulent flow with a reduced number of degrees of freedom. CVCS is also compared with a Fourier truncated simulation using a spectral cutoff filter where the number of retained Fourier modes is similar to the number of the wavelet coefficients retained by CVCS. It is shown that the wavelet representation is more suitable than the Fourier representation, especially concerning the probability density functions of vorticity and current density.


Introduction
Magnetohydrodynamic (MHD) turbulence is ubiquitous in astrophysical flows (Goldstein et al. 1995, Brandenburg andSubramanian 2005), as hydrodynamic (HD) turbulence is in our daily life. Turbulence has a large number of degrees of freedom, a wide range of dynamically active scales and strong nonlinearity. A key feature of turbulence is its strong spatial and temporal intermittency. The intermittency is attributed to coherent structures, e.g. vortex tubes for three-dimensional (3D) HD turbulence (e.g. Batchelor and Townsend 1949) and vorticity sheets and current sheets for 3D MHD turbulence (e.g. . Wavelet analysis is an efficient tool which yields a sparse multiscale representation of intermittent fields, because wavelets are well-localized functions in both space and scale. Thus, decomposing such fields into wavelet space keeps track of both their location and their scale. Wavelet techniques to analyze and simulate turbulent flows have been developed mainly for HD turbulence since the pioneering works (e.g. Farge and Sadourny 1989, Meneveau 1991, Yamada and Ohkitani 1991. Readers interested in the application of wavelets to turbulence may refer to review articles, e.g. Farge (1992), Schneider and Vasilyev (2010).
To reduce the degrees of freedom of turbulent flows and thus computational load in numerical simulations, a wavelet-based approach, called coherent vorticity simulation (CVS), was introduced for incompressible HD turbulence (Farge et al. 1999, Farge et al. 2003. CVS is based on the deterministic computation of the coherent flow evolution, using an adaptive wavelet basis while either modeling or neglecting the influence of the incoherent background flow. The underlying idea of CVS is to extract coherent vorticity out of HD turbulence at each time step. The coherent and incoherent vorticity fields are obtained by nonlinear filtering, i.e. thresholding the orthogonal wavelet coefficients of the vorticity field. The coherent vorticity is reconstructed from few wavelet coefficients whose modulus is larger than a given threshold, while the incoherent vorticity is obtained from the many remaining wavelet coefficients. The coherent and incoherent velocity fields are, respectively, computed from the coherent and incoherent vorticity fields, using the Biot-Savart kernel. The choice of the threshold is directly related to Donoho's (1993) criterion which supposes the incoherent flow to be Gaussian and decorrelated. At a given time instant, the coherent flow well preserves the statistics of the turbulent flow, while the discarded incoherent flow is like Gaussian and exhibits an energy equipartition. As the Reynolds number increases, the wavelet representation of 3D incompressible HD turbulence becomes more efficient in terms of the number of degrees of freedom .
To track the translation of the coherent vorticity and the generation of smaller scales, a safety zone is required. It corresponds to adding in wavelet space neighbors to the wavelet coefficients of the coherent vorticity at each time step. For 2D incompressible HD turbulence (Fro¨hlich andSchneider 1999, Schneider et al. 2006), 3D incompressible HD turbulent mixing layers (Schneider et al. 2005), and 3D incompressible HD homogeneous isotropic turbulence (Okamoto et al. 2011a), it was shown that CVS with the safety zone preserves the nonlinear dynamics and the statistical predictability, while in CVS without safety zone both are lost. Fully adaptive CVS will reduce not only the memory requirements, but also the CPU time of the computations of 3D incompressible HD turbulent flows, as shown for 3D compressible HD mixing layers (Roussel and Schneider 2010). The estimation of the numerical cost of fully adaptive CVS for the incompressible HD flows is discussed in (Okamoto et al. 2011a). The CVS approach is different from the filtering approach of large eddy simulation (LES), for the latter, see, e.g. Lesieur and Me´tais (1996), Meneveau and Katz (2000). In LES only the evolution of the large-scale flow is computed while modeling the influence of small-scale motion onto the large-scale motion.
Wavelet methods have also been applied for analyzing MHD turbulence. For example, Ishizawa and Hattori (1998) divided 2D MHD turbulent flows into two regions, using wavelet decomposition of the velocity and magnetic fields: one is the turbulent region showing Iroshnikov and Kraichnan spectra, i.e. the total energy spectra with the exponent of À3/2 (Iroshnikov 1964, Kraichnan 1965, and the other has spectra with À2 slope related to the existence of isolated current sheets. Recently, Yoshimatsu et al. (2011) examined the scale-dependent statistics of 3D MHD turbulence in the absence of a uniformly imposed magnetic field, and showed that the Lagrangian acceleration does not exhibit substantially stronger intermittency than the Eulerian acceleration. This behavior is in contrast to 3D HD turbulence where the Lagrangian acceleration shows much stronger intermittency than the Eulerian acceleration. Okamoto et al. (2011b) applied wavelet-based statistics to quasistatic MHD turbulence and quantified the flow anisotropy and its intermittency. By generalizing the coherent vorticity extraction method from HD turbulence, a method to extract coherent vorticity sheets and current sheets out of 3D MHD turbulence was introduced (Yoshimatsu et al. 2009). It was shown that only few degrees of freedom of both vorticity and current density, the coherent ones, well preserve the statistics of the total 3D MHD turbulence. Using this extraction method at each time step and neglecting or modeling the incoherent contributions, the time evolution of the coherent vorticity and current density can be simulated if a safety zone is introduced as in the case of HD turbulence. We call this simulation method coherent vorticity and current density simulation (CVCS).
The aim of this work is to develop a methodology of CVCS by generalizing CVS and then to show that the multiscale simulation method, CVS, has a great potential for application to other types of intermittent dynamics. To illustrate this, we perform CVCS for 3D forced homogeneous MHD turbulence in the absence of an imposed uniform magnetic field. Homogeneous turbulence is chosen here in order to demonstrate the efficiency of CVCS in the worst possible case where structures are spread all over physical space in contrast to inhomogeneous turbulence. Indeed, the wavelet representation is even more efficient for inhomogeneous flows, such as turbulent mixing layers, than for dealing with homogeneous flows. Forced turbulence is simulated in order to keep the small scales active, and to achieve the highest possible Reynolds number under the limitation of the available computational resources. To assess CVCS, the results are compared with direct numerical simulation (DNS) using the same maximal resolution. Comparing CVCS to a Fourier truncated simulation using a spectral cutoff filter, we examine whether the wavelet representation is more efficient with respect to the statistical predictability of turbulence than the Fourier representation. In this study we perform CVCS using pseudo-adaptive computations, as done for CVS in Schneider et al. (2005) and Okamoto et al. (2011a). The motivation is to get insight into the feasibility of fully adaptive CVCS computations, therefore we do not focus here on the CPU time.
The remainder of this article is organized as follows: in the next section, the CVCS method is presented after a short review of orthonormal wavelet representation. Section 3 shows the numerical method and results. Finally, conclusions are drawn in section 4.

Coherent vorticity and current density simulation method
We summarize the orthogonal wavelet decomposition of a generic 3D vector valued field. Then, we describe the coherent vorticity and current density simulation (CVCS) method, i.e. a method to simulate the evolution of coherent vorticity and current density, which is based on the wavelet filtered MHD equations.

Orthogonal wavelet decomposition of 3D vector field
We consider a 3D 2-periodic vector field v( Here, J is the number of octaves in each space direction of the Cartesian coordinates, The vector field v, having a mean value " v (which vanishes in the flows studied here), can be decomposed into an orthogonal wavelet series: where the index set of the wavelet coefficients L is . . , 7, j ¼ 0, . . . , J À 1, and i n ¼ 0, . . . , 2 j À 1 ðn ¼ 1, 2, 3Þ È É : The multi-index denotes the scale 2 Àj and the position 2 Á 2 Àj i ¼ 2 Á 2 Àj (i 1 , i 2 , i 3 ), and denotes the seven directions of the wavelets. The index set L can be seen as the octree representation of the orthogonal wavelet coefficients. We have 7 Â 2 3j wavelet coefficients at scale 2 Àj for each component of v. A family of 3D wavelets , (x) is generated by dilation and translation of the 3D mother wavelet (x) based on a tensor product construction: Here (x) is the one-dimensional mother wavelet, and (x) is the one-dimensional scaling function. The family of , yields an orthogonal basis of not only L 2 (R 3 ) but also of L 2 (O) through the application of a periodization technique (see, e.g. Mallat 2009).
Note that , (x) is well-localized in space x 2 R 3 , oscillating, smooth, and the spatial average of , (x) vanishes for each index. Owing to orthogonality, the wavelet coefficients e v , are given by e v , ¼ v, , , where hÁ, Ái denotes the L 2 -inner product, defined by h f, gi ¼ The coefficients measure the fluctuations of v at scale 2 Àj and around position 2 Á 2 Àj i for each of the seven directions ¼ 1, . . . , 7. The N À 1 wavelet coefficients e v , and one mean value " v are efficiently computed from the N ¼ 2 3J grid point values of v using the fast wavelet transform which has linear computational complexity. In this work, the compactly supported Coiflet wavelets with filter width 12 are used. For details on the orthogonal wavelet transform, we refer readers to, e.g. Farge (1992) and Mallat (2009).

Coherent vorticity and current density extraction method
We present a summary of the method to extract both coherent vorticity sheets and coherent current sheets out of 3D MHD turbulence, which was introduced by Yoshimatsu et al. (2009). This method is based on orthogonal wavelet decompositions of vorticity x ¼ ; Â u and current density j ¼ ; Â b. Here u is the velocity field and b is the magnetic field. The latter is normalized by ( 0 0 ) 1/2 , where 0 is the permeability of free space and 0 is the fluid density.
Thresholding the wavelet coefficients e x , and e j , at a given time instant, with a thresholding function we separate the coefficients into coherent and incoherent ones, i.e. e v c ¼ T ðe v , Þ and e v i ¼ e v , Àe v c , respectively. The value of Tð e x , Þ can be different from that of Tð e j , Þ.
The coherent vorticity, x c , the incoherent vorticity, x i , the coherent current density, j c , and the incoherent current density, j i , are then reconstructed by the inverse wavelet transform. These fields satisfy Kinetic and magnetic enstrophies are conserved, i.e.
The choice of the thresholds is motivated by denoising theory (Donoho 1993, Donoho andJohnstone 1994). The values of the thresholds are given by Tð e x , Þ ¼ fð4=3ÞZ u i ln Ng 1=2 and Tð e j , Þ ¼ fð4=3ÞZ b i ln Ng 1=2 , which depend only on the incoherent enstrophies Z u i and Z b i (which are a priori unknown) and the maximal resolution N. Parseval's identity implies that the enstrophy can be computed directly in wavelet space. We use the threshold values, T 0 ð e x , Þ and T 0 ð e j , Þ, from the total enstrophies, Z u and Z b , as the first estimates of Z u i and Z b i , respectively, i.e. T 0 ð e x , Þ ¼ fð4=3ÞZ u ln Ng 1=2 and T 0 ð e j , Þ ¼ fð4=3ÞZ b ln Ng 1=2 . Then we split x and j into coherent and incoherent fields by thresholding their wavelet coefficients with T 0 ð e x , Þ and T 0 ð e j , Þ, respectively. The obtained incoherent enstrophies provide new estimates of Z u i and Z b i , from which new threshold values T 1 ð e x , Þ and T 1 ð e j , Þ are computed. The above procedure could then be iterated until convergence. The coherent vorticity and current density extracted by the wavelet thresholding with T 1 ð e x , Þ and T 1 ð e j , Þ are sufficient to preserve the statistics of the total field at a given time instant and yield good compression. In contrast, the coherent vorticity and current density extracted by the wavelet filtering with T 0 ð e x , Þð4T 1 ð e x , ÞÞ and T 0 ð e j , Þð 4 T 1 ð e j , ÞÞ do not sufficiently preserve the statistics of the total field.
Biot-Savart's relations, u ¼ À ; Â (r À2 x) and b ¼ À ; Â (r À2 j), are used to obtain the corresponding induced velocity and magnetic fields, respectively. By construction the total fields are obtained by summation, For the kinetic and magnetic energies, we have The incoherent velocity and magnetic fields are like Gaussian and exhibit an energy equipartition.

Wavelet filtered MHD equations
In order to follow the time evolution of MHD turbulence, we should apply the method to extract coherent vorticity and current density at each time step. However, to track the evolution, adding a safety zone in the vicinity of the coherent wavelet coefficients of vorticity and current density, e x c and e j c , is inevitable, as it is the case for CVS. Thus, the CVCS method should retain not only e x c and e j c but also their neighbor wavelet coefficients at each time step.
Let L n c,x and L n c, j be the index sets of the coherent wavelet coefficients of x c and j c at a given time instant t n , respectively. We then consider the union of these sets denoted by L n cþ , i.e. L n cþ ¼ L n c,x [ L n c, j . Subsequently, the index set L n cþ is expanded by adding neighbors in direction, space and scale. The zone consisting of the neighbors is called the safety zone. There are 2 or 3 in direction, 26 adjacent neighbors in space, and 8 in scale for a given wavelet coefficient, as illustrated in figure 1. The expanded index set, denoted by L n Ã , allows to track not only the translation of x c and j c but also the production of finer scales due to their nonlinear interactions. To take advantage of the safety zone as much as possible, we keep the incoherent wavelet coefficients belonging to the safety zone, i.e. the difference set L n Ã nL n cþ , at any time. We neglect the remaining large majority of the wavelet coefficients, which are not elements of L n Ã . One should verify whether the width of the safety zone, which depends on the time increment, is sufficient to preserve the flow dynamics well, because the incompressibility conditions given by equations (8) and (9) imply that local changes in the flow could become global. For HD turbulence, the statistics of the reference DNS are shown to be well preserved by CVS if a safety zone is provided (Schneider et al. 2005, 2006, Okamoto et al. 2011a. These studies suggest that information which travels further than the adjacent wavelet coefficients are not significant and can thus be neglected.
The coherent fields including the safety zone, x cÃ and j cÃ , are obtained from the wavelet coefficients belonging to L Ã . The remaining vorticity and current density, x iÃ and j iÃ , are given by x iÃ ¼ x À x cÃ and j iÃ ¼ j À j cÃ , respectively. The corresponding velocity and magnetic fields, u cÃ , u iÃ , b cÃ and b iÃ , which are induced by x cÃ , x iÃ , j cÃ and j iÃ , are obtained from the Biot-Savart relation. Although x cÃ and j cÃ are not perfectly solenoidal in CVCS, the solenoidal conditions are ensured by taking the curl of u cÃ and b cÃ . In section 3.2, we will confirm that the discrepancies between x cÃ and ; Â u cÃ and between j cÃ and ; Â b cÃ are indeed negligible. Substituting into the Navier-Stokes equations and the induction equations, and then neglecting the terms including x iÃ , u iÃ , j iÃ and b iÃ , we obtain the evolution equations for u cÃ and b cÃ : ; · u cÃ ¼ 0, where f u and f b are external solenoidal forces imposed on large-scale flow, P is the pressure, is the kinematic viscosity, and is the magnetic diffusivity. The evolution of x cÃ and j cÃ can be obtained by taking the curl of equations (6) and (7). A flowchart of the CVCS algorithm is shown in figure 2.

Numerical method
We performed five numerical computations of forced incompressible MHD turbulence without mean magnetic field in a 2 periodic box; one DNS computation as reference, three CVCS computations (CVCS0, CVCS1 and CVCS2) and one Fourier truncated simulation (FT0) using a spectral cutoff filter. In FT0, only Fourier coefficients of the velocity and magnetic fields for jkj k c are retained, while coefficients with wavenumbers larger than k c are set to zero. Here, k ¼ jkj and k c is the cutoff wavenumber of the filter where k is the wave vector. The magnetic Prandtl number Pr is set to 1, i.e. ¼ . CVCS0 performs wavelet thresholding with T 0 ð e x , Þ and T 0 ð e j , Þ (without iteration), while CVCS1 and CVCS2 use threshold values T 1 ð e x , Þ and T 1 ð e j , Þ applying one iteration. The thresholds for e x , and e j , depend on the kinetic and magnetic enstrophies, respectively, and thus the threshold values vary in time. CVCS0 and CVCS1 include the safety zones, while CVCS2 has no safety zone. An overview on the different CVCS computations is given in table 1. We perform these CVCS using pseudo-adaptive computations, as mentioned in section 1. We do not employ any subgrid scale model in FT0, because the influence of the incoherent flow onto the coherent flow is neglected for the present CVCS either. The cutoff wavenumber k c of the spectral cutoff filter is chosen such that the percentage of the retained Fourier coefficients is almost the same as for the retained wavelet coefficients in CVCS0. The numerical code uses a classical Fourier pseudo-spectral method based on the Elsa¨sser variables formulation with a fourth-order Runge-Kutta scheme for time marching. The variables z AE are defined by z AE ¼ u AE b. In each time step, we compute the coherent vorticity and current density wavelet coefficients, e x c and e j c , by applying the coherent vorticity and current density extraction method, add the corresponding safety zone, and invert the curl operator to compute the coherent velocity and magnetic field, u cÃ and b cÃ . Thus we obtain the corresponding Elsa¨sser variables z AE cÃ ¼ u cÃ AE b cÃ , which are advanced in time. The aliasing errors are removed by a phase-shift method, which keeps all the Fourier modes satisfying k < k max , where k max ¼ 2 1/2 N 1/3 /3. Solenoidal random forces introduced by Yoshida and Arimitsu (2007) are imposed only at large scale, i.e. in the wavenumber range 1 k < 2.5. The resolution is N ¼ 256 3 , ¼ ¼ 10.2 Â 10 À4 , and the time increment Át ¼ 4.2 Â 10 À3 . For both forces, f u and f b , the correlation time and the intensity are, respectively, set to 1.8 and 0.9 Â 10 À3 . Note that they have the same time history in all the presented computations.
The computations are integrated over the time period 3.31, where ¼ L/u 0 , u 2 0 ¼ 2E u =3 and L is the integral length scale defined by L ¼ =ð2u 2 0 Þ R k max 0 E u ðkÞ=k dk and k max is the maximum wavenumber. The kinetic energy spectrum E u (k) is defined by E u (k) ¼ P kÀ1/2 jqj<kþ1/2 jF [u](q)j 2 /2. Here, F [Á] is the Fourier transform of Á. Typical physical quantities at the final time t f ¼ 3.31 are summarized in table 2.
Notes: The symbols * and Â denote the computations with and without the safety zone, respectively. The threshold T 1 stands for a threshold with one iteration, and T 0 for a threshold without iteration. Note that T 0 > T 1 . We use the same initial flow field in all computations. This field was obtained by a preceding DNS computation of 3D incompressible homogeneous MHD turbulence with the kinetic Taylor micro-scale Reynolds number R u ¼ 86, the magnetic Taylor micro-scale Reynolds number R b ¼ 124 and k max IK ¼ 1.76. The total energy dissipation rate hi is statistically quasi-stationary. Here, R u ¼ u 0 u = and R b ¼ b 0 b =, the Iroshnikov and Kraichnan microscale IK is defined by ( 2 b 0 /hi) 1/3 , where the kinetic Taylor microscale u ¼ ð15u 2 0 =h u iÞ 1=2 , the magnetic Taylor microscale b ¼ ð15b 2 0 =h b iÞ 1=2 , b 0 ¼ (2 " E b /3) 1/2 , h u i and h b i are the kinetic and magnetic energy dissipation rates, respectively. The initial velocity and magnetic fields of the preceding DNS, which are randomly generated, have the following constraints; E u (k) ¼ E b (k) / k 2 exp(À2k 2 /25), the total cross helicity and magnetic helicity are nearly zero, and both the initial total kinetic and magnetic energies are 0.5.

Assessment of CVCS
In this subsection, we study the feasibility of CVCS and assess its quality in comparison to DNS and the Fourier truncated simulation.
We first describe the percentage C of the wavelet coefficients retained by CVCS, which is defined by C ¼ 100N Ã /N. Here N Ã is the number of the wavelet coefficients which are elements of L Ã , in other words, the number of wavelet coefficients retained by the coherent vorticity and the coherent current density including the safety zone. Figure 3 shows C versus normalized time t/. In each CVCS, C is almost independent of time after a transient decay at early times, for t/ 0 0.1. CVCS1 retains most coefficients (C $ 33[%]), while CVCS0 exhibits a smaller percentage, C $ 13[%]. This is consistent with the fact that CVCS0 employs larger threshold values for the wavelet coefficients of vorticity and current density than CVCS1. For CVCS2, which has no safety zone, the percentage C rapidly decreases to much smaller value (about 0.9%). For FT0 we have chosen a cutoff wavenumber k c ¼ 80 such that the percentage of retained Fourier coefficients yields about 13%, which is thus slightly larger than the percentage of CVCS1 for t/ > 1.   This observation is in contrast to what is found in CVS of 3D homogeneous incompressible HD turbulence (Okamoto et al. 2011a). The CVS computations show a statistically similar picture of entangled vortex tubes as in DNS. However, the position of these intense vorticity regions in CVS is completely different from those in DNS because of the flow sensitivity. Considering finally the CVCS2 computation, where no safety zone has been added, we find that the sheet-like structures are not preserved and thus it is suggested that the flow dynamics is lost. It was confirmed that impressions obtained from the visualization are the same as those at different time instants t t f , say, t ¼ t f /2 (figure omitted).
Let us define the ratio of a quantity ''X'' in CVCS to that of DNS by R(X), in order to quantify how much kinetic energy, E u , magnetic energy, E b , kinetic enstrophy, Z u , and magnetic enstrophy, Z b , are retained by CVCS compared to DNS. Figures 6(a)-(d) illustrate the time evolution of 100R(E u ), 100R(E b ), 100R(Z u ) and 100R(Z b ), respectively. In figures 6(a) and (b), it is seen that 100R(E u ) and 100R(E b ) are almost 100% for CVCS0 and CVCS1, i.e. the kinetic and magnetic energies are in excellent agreement with those in DNS, during the time evolution period studied here. The differences of the energies among these CVCS and DNS are small and less than 1%. In contrast, the evolution of the ratios in CVCS2 differs by 3 % for kinetic energy and by 6% for magnetic energy. FT0 yields reasonable predictability for the time evolution of the energies with values being slightly above the ones of DNS.
In figures 6(c) and (d), it is seen that CVCS1 yields large values of enstrophy for both Z u and Z b than CVCS0, because CVCS1 has more retained wavelet coefficients  Figure 7 shows the flatness versus t/. Both of the flatness values of F ! and F j for CVCS0 and CVCS1 are a little larger than the corresponding values of DNS. In contrast, the values in CVCS2 are much smaller than F ! and F j in DNS. The flatness values F ! and F j for FT0 are also much smaller than those for DNS. Therefore, we find that the wavelet representation is more suitable for predicting the fourth-order turbulence statistics than the Fourier representation.
The probability density functions (PDFs) of velocity and magnetic fields for all computations agree well with each other (figures 8(a) and (b)). All PDFs exhibit a Gaussian shape. However, this is not the case for the PDFs of vorticity and current density, as shown in figures 8(c) and (d). CVCS0 and CVCS1 well preserve the PDFs of the vorticity and current density obtained by DNS, whose tails have stretched exponential shapes. On the other hand, in the PDFs of the vorticity and current density for CVCS2, their supports are significantly narrower compared to those for DNS, respectively. The supports of the PDFs of vorticity and current density for FT0 are also reduced compared to DNS and CVCS0, although the supports for FT0 are wider than for CVCS2.
The kinetic and magnetic energy spectra, E u (k) and E b (k), are plotted in figure 9, to assess the quality of CVCS in terms of two-point statistics. We find that, for CVCS0 and CVCS1, both of the energy spectra agree well with those for DNS, respectively. CVCS0 and CVCS1 preserve the DNS spectra in the range k IK 0 0.5, while for larger wavenumbers some departure, which is less pronounced for CVCS1 than for CVCS0, can be observed. For CVCS2 we observe a strong departure for k IK 0 0.2. In contrast to CVCS, FT0 tends to pile up energy in the spectra near the cutoff wavenumber k c .
Here, x ' and j ' are '-th components of x and j, respectively.
A consequence of this effect is that the kinetic and magnetic enstrophies in figure 6 have values of about 100%.
The total energy flux Å (k) in figure 10 confirms that the nonlinear dynamics are well retained by CVCS0, CVCS1 and FT0. Here Å(k) is defined by (a) (b) Figure 9. (color online) (a) Kinetic and (b) magnetic energy spectra E u (k) and E b (k) at t ¼ t f . The wavenumber k is normalized by IK of the DNS at t ¼ t f .
Again we see that CVCS2 exhibits substantial departure in terms of the energy flux.
To further quantify the difference of the simulations, we consider the spectra of the kinetic and magnetic difference fields, u À u DNS and b À b DNS ( ¼ CVCS0, CVCS1, FT0). Here u and b express the velocity and magnetic field of the computation , respectively, u DNS the velocity for DNS, and b DNS the magnetic field for DNS. The spectra of the kinetic and magnetic difference fields are respectively defined as In figure 11, we see that D u CVCS0 ðkÞ and D b CVCS0 ðkÞ are comparable to the kinetic and magnetic energy spectra E u (k) and E b (k) for DNS in the range k IK 0 0.5, respectively. Around k IK $ 0.2, these spectra of the difference fields take their maximum: the values are about 10% of E u (k) and E b (k) for DNS. For the difference spectra of CVCS1, we have D u CVCS1 ðkÞ 5 D u CVCS0 ðkÞ and D b CVCS1 ðkÞ 5 D b CVCS0 ðkÞ, as expected. Concerning FT0, we find D u FT0 ðkÞ 5 D u CVCS0 ðkÞ and D b FT0 ðkÞ 5 D b CVCS0 ðkÞ for k IK 9 1, because the Fourier modes of the difference fields increase for k k c more rapidly in CVCS0 than in FT0. CVCS0 removes some contributions in the range k k c by filtering out the incoherent contribution, while CVCS0 retains some contributions in the coherent flow for k > k c . The difference spectra D u and D b grow in time at each wavenumber until they become comparable to E u (k) and E b (k) at least up to t ¼ t f , respectively (figure omitted), as is the case studied by Yamazaki et al. (2002) who examined the effects of wavenumber truncation in DNS of 3D homogeneous incompressible HD turbulence in a periodic box using the Fourier spectral method. Thus, we confirm that CVCS cannot preserve the deterministic predictability, as any truncated method, although we can see from the visualization in figures 4 and 5 that the positions of the vorticity sheets and current sheets are well preserved by CVCS0 and CVCS1 during the present simulation time period. It is suggested that the positions are not retained at later times. Furthermore, from figures 4, 5, and 11 we conjecture that the dynamics in the range k IK 9 0.2 play a key role in determining the position of the sheet-like intense structures.
Finally, let us examine the compressible contributions of the coherent vorticity x cÃ and coherent current density j cÃ and confirm that they are not crucial in CVCS with the safety zone. Figure 12 shows the spectra of x cÃ and j cÃ , their divergence-free parts, x 0 cÃ and j 0 cÃ , and their compressible contributions, ; ! and ; j at the final time in the case of CVCS0. Here x 0 cÃ ¼ x cÃ þ ; ! and j 0 cÃ ¼ j cÃ þ ; j . We observe that the spectra of x cÃ and j cÃ agree well with those of x 0 cÃ and j 0 cÃ , respectively. The compressible contributions are observed mainly in the dissipative range. However, they are some orders of magnitude below the corresponding incompressible parts. The relative errors of the compressible contributions are estimated by 100 Â hj; ! j 2 i/hjx cÃ j 2 i $ 0.17 [%], and 100 Â hj; j j 2 i/hj j cÃ j 2 i $ 0.16 [%]. Therefore, the compressible contributions are (a) (b) Figure 11. (color online) Spectra of (a) kinetic and (b) magnetic difference fields at t ¼ t f : the dashed curves denote D u CVCS0 and D b CVCS0 , the dashed dotted curves D u CVCS1 and D b CVCS1 , and the dotted ones D u FT0 and D b FT0 . As references, the kinetic and magnetic energy spectra for DNS at t ¼ t f are plotted as thick solid curves, respectively. The wavenumber k is normalized by IK of the DNS at t ¼ t f .
(a) (b) Figure 12. (color online) Spectra of the coherent field a cÃ , its incompressible contribution a 0 cÃ , and its compressible one ; at t ¼ t f for CVCS0: (a) for the coherent vorticity, i.e. a ¼ x, and (b) for the coherent current density, i.e. a ¼ j. The wavenumber k is normalized by IK of the DNS at t ¼ t f . negligible in the CVCS. Note that CVCS1 has smaller relative errors than CVCS0, because the former retains more wavelet coefficients than the latter (figure omitted).

Conclusions
We have developed the CVCS method to track the time evolution of coherent vorticity and current density for 3D incompressible MHD turbulence, generalizing the CVS method for 3D incompressible HD turbulence previously introduced by . CVCS is based on the orthogonal wavelet decomposition of the vorticity and current density fields at each time step. Each of the coherent fields corresponds to few wavelet coefficients whose modulus is larger than a threshold, and the incoherent fields are obtained from the remaining large majority of the coefficients. In order to obtain the flow evolution of these coherent fields, one needs to add a safety zone in wavelet space. The wavelet coefficients which are computed in the next time step consist of the coherent ones and their neighbors, i.e. the safety zone. The remaining incoherent wavelet coefficients are discarded in the present CVCS computations. The presented results show that discarding these coefficients is sufficient to model turbulent dissipation.
CVCS of 3D forced homogeneous incompressible MHD turbulence in the absence of an imposed uniform magnetic field has been performed to assess its quality. CVCS was compared with a reference DNS at an initial kinetic Taylor microscale Reynolds number R u ¼ 86, and an initial magnetic one R b ¼ 124. These computations were performed at resolution N ¼ 256 3 . Homogeneous turbulence was chosen to demonstrate the feasibility and efficiency of CVCS in the most challenging case.
We found that the statistics of the reference DNS are well preserved by CVCS if the safety zone is provided. CVCS1, which includes the safety zone, presents the best statistical predictability, but the number of the degrees of freedom retained by CVCS1 is only reduced by a factor three in comparison to DNS. To improve the computational efficiency, we performed CVCS0, which uses larger threshold values than CVCS1. CVCS0 has the same type of safety zone as CVCS1. However, the former reduces the number of degrees of freedom by a factor eight, in comparison to DNS. Thus, it was verified that information which travels further than the adjacent wavelet coefficients is not crucial to track the evolution of the nonlinear dynamics, since the CVCS computations with safety zone well retain the statistics of DNS. In contrast, CVCS2 without any safety zone loses both the flow dynamics and the flow statistics. The small divergence contribution in the coherent vorticity field was confirmed to be negligible. Therefore, it is demonstrated that the CVS method has also a great potential to be applied to simulations of flows whose intermittency is different from HD turbulence.
As any truncated method the deterministic predictability of turbulent flows cannot be preserved by CVCS owing to the flow sensitivity. Flow visualization suggests that CVCS better preserves the positions of intense vorticity sheets and current sheets during the simulation, compared to CVS of 3D incompressible homogeneous HD turbulence in which the positions of vortex tubes are completely different from those for DNS.
We have also compared CVCS with a Fourier truncated simulation. By the use of the spectral cutoff filter, the number of the modes retained by the Fourier truncated simulation is set to be similar to the number of wavelet coefficients retained by CVCS0.
For the sake of comparison no subgrid scale model was added to the Fourier truncated simulation, because the influence of the incoherent flow onto the coherent flow is not modeled in CVCS either. We found that the wavelet representation better predicts the turbulence statistics such as the vorticity and current density PDFs, than the Fourier representation.
In future work, the development of fully adaptive CVCS computations, especially for MHD turbulence at much higher Reynolds number will be pursued. For CVCS of inhomogeneous flows in complex geometries, a hybrid method between the volume penalization approach for MHD flows, proposed (Neffaa et al. 2008, and CVCS is promising. A generalization of the CVCS methodology to anisotropic MHD turbulence under a substantial uniformly external magnetic field with and without rotation, density stratification, the Hall effect, etc., would also be interesting. The safety zone is related to physical processes such as sweeping, energy cascade, etc. Indeed, the CVCS approach is based on the assumption of locality of energy transfer in wavelet space, which allows an extrapolation of the flow dynamics in time. Adding a safety zone, i.e. including adjacent wavelet coefficients, allows to capture the translation of coherent vorticity and current sheets in space and direction, and the generation of finer scales due to nonlinear interaction. For substantially anisotropic turbulence, a significant reduction of the width of the safety zone is expected, because the energy cascade and the small scale motion are suppressed in the direction of the external magnetic field. In some situations, dynamos may occur, and then kinetic and magnetic helicities may play significant roles in the dynamics and statistics. It is conjectured that the above arguments on the safety zone still hold and that the dynamics and statistics are well preserved by the coherent flow of the MHD turbulence. In Farge et al. (1999), it was shown that CVS tracks well the coherent flow evolution in the regime of an inverse energy cascade for 2D HD turbulence. However, examination of this conjecture for MHD is beyond the scope of this article.
for hospitality during the 2010 CEMRACS summer program on 'Numerical modeling of fusion'.