A WELL-POSED KELVIN–HELMHOLTZ INSTABILITY TEST AND COMPARISON

Kelvin–Helmholtz instability (KHI) is the name given to the primary instability that occurs when velocity shear is present within a continuous fluid or across fluid boundaries. The shear is converted into vorticity that, subject to secondary instabilities, cascades generating turbulence. The KHI is one of the most important hydrodynamical instabilities and plays a significant role in various parts of astrophysics. It is believed to be responsible for additional mixing in differentially rotating stellar interiors (Brüggen & Hillebrandt 2001), and to keep a finite-thickness layer of dust around the midplane of protoplanetary disks (Dubrulle et al. 1995; Johansen et al. 2006). It also contributes to convective mixing in any deep stellar interior at the stiff convective boundaries, for instance, in asymptotic giant stars (Herwig 2006) or novae (Casanova et al. 2011). Moreover, KHI can lead to the destruction of cool gravitationally bound objects moving in a hot ambient medium (Murray et al. 1993) such as galaxies in the intracluster medium (Nulsen 1982; Mori & Burkert 2000), substellar companions engulfed by a giant star (J.-C. Passy et al., in preparation) and comets entering a planetary atmosphere (Mac Low & Zahnle 1994). KHI plays a role in the interactions of the magnetopause and solar wind (Miura & Pritchett 1982) and has been observed in the solar corona (Ofman & Thompson 2011). In order to understand these phenomena and their implications, it is therefore important to define a well-posed method to quantify how accurately KHI can be modeled by different numerical techniques.

Verifying the correct treatment of KHI has attracted increased interest following the conclusions made by Agertz et al. (2007) including vigorous discussions of KHI in Lagrangian schemes. The main conclusion reached was that smoothed particle hydrodynamics (SPH) fails to resolve KHI due to a surface tension effect between the SPH particles at the shear interface. However, the test was done at a sharp shear and contact discontinuity. Price (2008) attempted to address the problem with KHI growth from a sharp contact discontinuity in Agertz et al. (2007) by adding an artificial thermal conductivity to SPH. A prescription achieving a similar end by adding a diffusion motivated by a subgrid turbulence model to SPH was proposed by Wadsley et al. (2008). In a case where traditional SPH largely fails to reproduce KHI at a sharp interface, Cha et al. (2010) demonstrated that using a Godunov–SPH formulation with zeroth- and first-order consistency that growth of KHI can be obtained. Using a Voronoi-mesh-based scheme Heß & Springel (2010) showed improvement over SPH in a sharp contact discontinuity KHI test, but compared the compressible solution to the growth rate for an incompressible flow and did not perform a convergence study. With the AREPO Voronoi-mesh Godunov code Springel (2010a) ran a sharp contact discontinuity KHI test and pointed out the difference seen in the secondary instabilities developed when the mesh was given a motion following the flow (a quasi-Lagrangian motion) or kept fixed. In Springel (2011) the same code is used for a extended discussion, with both a sharp contact discontinuity KHI test and a smooth transition test, but comparing both of the compressible results to the growth rate for an incompressible sharp contact continuity initial test. Read et al. (2010) pointed out the zeroth-order inconsistency in SPH, and designed a kernel to minimize these effects, achieving better qualitative results on a sharp contact discontinuity KHI test. Zeroth-order inconsistency is the inability of SPH interpolation to reproduce a constant function at any finite resolution (Liu et al. 1995; Dilts 1999; Liu et al. 2003; Fries & Matthies 2004; Quinlan et al. 2006). A quantitative analysis of the growth of KHI from a sharp contact discontinuity was performed by Junk et al. (2010) with SPH and grid-based Godunov codes. With a focus on SPH, Valcke et al. (2010) argued that a sharp contact discontinuity KHI test is not ideal and propose an alternative SPH smoothing kernel which yields improved results. Hubber et al. (2011) compared qualitatively a grid-based method and SPH using both a cubic and quintic kernel on a sharp contact discontinuity KHI test, finding that the choice of a quintic kernel improved the SPH results significantly. One of the only well-posed convergence tests for KHI was done in Robertson et al. (2010), but that was in a study of Galilean invariance restricted to fixed-mesh schemes. The test by Springel (2010b) is a well-posed problem influenced by Robertson et al. (2010), but the evaluation of the SPH result was done by comparison to an analytic solution for a sharp transition initial condition and incompressible flow in an infinite domain, not for the problem posed with a softened transition in compressible flow in a finite periodic domain.

The commonly used solution for the KHI growth rate in numerical tests is for a sharp transition at the shear interface (Chandrasekhar 1961, Section 101). However, for numerical approximations the interface should be smoothed to yield an initial value problem with finite spatial derivatives, as argued by Robertson et al. (2010). For a sharp interface, the initial approximation of the derivatives across the interface does not converge with resolution. To obtain convergence, a smooth interface must be used. We pose the problem in such a way that the analytic result for the incompressible limit is known for an infinite domain, as this type of analytic result is usually used to compare numerical results. However, a difficulty with Kelvin–Helmholtz problems is that the unstable modes are global, so solutions in a finite periodic domain, as commonly tested, are different from than the solution in an infinite domain. To circumvent this difficulty, we perform an exhaustive convergence study to establish a fully compressible nonlinear solution with a very small and rigorously derived uncertainty.

In the following discussion, we will refer to different types of discretizations used for numerical solutions to the chosen governing equations of hydrodynamics or magnetohydrodynamics. To clarify, most fixed grid (or structured mesh) codes use a square Eulerian grid as a basis for either a point value (values of the fields at grid points) or volume-average (average volume of the field in a grid cell) discretization. The distinction here is that a finite-volume scheme can be arranged to solve the integral form of the governing equations, and the point-value discretization can only solve the differential form of the equations. Unstructured mesh discretizations do not pose the restriction of the discretization mesh being a regular grid; however, the nodes are logically connected by edges, and the mesh cells form a tessellation on the computational volume. Moving-mesh Voronoi tessellation discretizations have begun to appear in astrophysical applications (Springel 2010a; Duffell & MacFadyen 2011). These are a case of an unstructured mesh, where the mesh is defined by the Voronoi tessellation of a set of mesh generating points. The mesh cells may be used to define a finite-volume discretization. When the mesh generating points are allowed to move in time to make a moving-mesh discretization, the Voronoi mesh will be recalculated on the new point distribution at every step yielding a new set of cells and new mesh edges connecting them. The mesh movement can be arbitrary, but a quasi-Lagrangian mesh movement is a particularly good choice as this minimizes numerical errors associated with advection across the grid. It is also possible to define meshless discretizations that represent the fields on a set of points or particles without specifying a set of mesh edges to connect these points. Meshless discretizations then do not form a strict tessellation of the computational volume. These discretizations are commonly used to define Lagrangian methods, in the sense that the points or particles are comoving with the fluid. If finite-mass particles are used, then the particles form a partition of the total computational mass, and the form of discretization used in SPH is obtained and the integral form of the governing equations can be solved. However, if the meshless points carry only field values at those points, then a point-value-type discretization is obtained and the differential form of the governing equations must be solved.

In Section 2, we give the problem setup used, and in Section 3 we discuss the methods used to extract measured quantities from the results. The codes used in this paper are listed in Section 4. The detailed convergence study used to generate the reference compressible solution is presented in Section 5. The results and comparison from several codes with different underlying algorithms and discretizations are presented in Section 6. We discuss the various results and implications in Section 7. Extended discussion of the SPH results and analysis of extra experiments are presented in Section 8. In Section 9, we discuss secondary instabilities arising from the problem setup in this work, and the difficulty of determining if they are produced in a physically meaningful manner. Our conclusions are summarized in Section 10.

Our motivation in choosing the initial condition is that the initial conditions are smooth, reflect as closely as possible a configuration that can be treated analytically, and can be represented easily in a wide variety of codes. In all codes, we will solve the inviscid compressible Euler equations. The setup we use is chosen to be a periodic version of that used in the analysis of KHI in Wang et al. (2010): the domain is 1 unit by 1 unit in the x- and y-directions if two-dimensional, and an arbitrary thickness in the z-direction if needed for a three-dimensional code. Runs with resolutions of 128 × 128, 256 × 256, and 512 × 512 cells, or equivalent were used in the comparison. All boundaries are periodic. The initial condition is smooth and periodic, as illustrated in Figure 1. The density is given by

$$\begin{equation} \rho = \left\lbrace \begin{array}{@{} ll}\rho _1 - \rho _m e^\frac{y-1/4}{L} &\mbox{ if } 1/4>y\ge 0\\ \ms \rho _2 +\rho _m e^\frac{-y+1/4}{L} &\mbox{ if } 1/2>y\ge 1/4\\ \ms \rho _2 +\rho _m e^\frac{-(3/4-y)}{L} &\mbox{ if } 3/4>y\ge 1/2\\ \ms \rho _1 -\rho _m e^\frac{-(y-3/4)}{L} &\mbox{ if } 1>y\ge 3/4 \end{array} \right. \end{equation} \tag{ 1 }$$

where

$$\begin{equation} \rho _m = (\rho _1 - \rho _2)/2 \end{equation} \tag{ 2 }$$

with ρ₁ = 1.0, ρ₂ = 2.0, and the smoothing parameter L = 0.025. The x-direction velocity is given by

$$\begin{equation} V_x = \left\lbrace \begin{array}{@{} ll}U_1 - U_m e^\frac{y-1/4}{L} &\mbox{ if } 1/4>y\ge 0\\ \ms U_2 + U_m e^\frac{-y+1/4}{L} &\mbox{ if } 1/2>y\ge 1/4\\ \ms U_2 + U_m e^\frac{-(3/4-y)}{L} &\mbox{ if } 3/4>y\ge 1/2\\ \ms U_1 - U_m e^\frac{-(y-3/4)}{L} &\mbox{ if } 1>y\ge 3/4\\ \end{array} \right. \end{equation} \tag{ 3 }$$

where

$$\begin{equation} U_m = (U_1 - U_2)/2 \end{equation} \tag{ 4 }$$

with U₁ = 0.5, U₂ = −0.5, and L as in the density so that the smooth transition in density and velocity occurs over the same interval. The background shear is perturbed by adding some velocity in the y-direction with the form

$$\begin{eqnarray} &&V_y = 0.01\sin (4\pi x). \end{eqnarray} \tag{ 5 }$$

An ideal gas equation of state with γ = 5/3 is used. The internal energy is set such that pressure is initially uniform with value 2.5. The problem is run till at least time t = 1.5. Analysis is done on snapshots spaced at a minimum of Δt = 0.02. However, in most cases, the snapshots will not be spaced exactly as codes often do output or analysis on an approximate interval, e.g., at the first time step after the specified snapshot or analysis time. The test can be run in two dimensions in a structured grid code, but for unstructured meshes or mesh-free methods two-dimensional and three-dimensional simulations may yield slightly different results depending on how the resolution elements are arranged in the initial condition. For unstructured mesh methods and meshless methods the results will differ for a disordered node distribution and a regularly gridded one.

Zoom In Zoom Out Reset image size

To quantitatively describe the growth of the KHI, we use two measurements, the amplitude of the y-velocity mode of the instability, and the maximum y-direction kinetic energy density. These two quantities are a useful pair, as the mode amplitude is a smoothed quantity and the maximum y-direction kinetic energy density is very sensitive to noise in the computed velocity field.

As a loose guide, the analysis of Wang et al. (2010) treats a non-periodic incompressible version of the problem studied in this paper. Their linear perturbation theory yields growth rates for the two quantities studied in this work, in the infinite domain and incompressible flow limit. However, as we run our test in periodic boundaries with a compressible flow we must go further than their analysis.

The maximum y-direction kinetic energy is the simplest of the two quantities to compute. This quantity is the maximum value of 1/2ρV²_y computed for all resolution elements (cells, points, or particles) in the computation volume at each time. In the non-periodic, incompressible limit, the growth of this quantity should be ∝exp (2 × 4.384 × t) (Wang et al. 2010; Equation (18)). In practice, the growth will start from a finite perturbation, will reflect erroneous velocities occurring both at the interface due to unbalanced pressures at the cell scale, and any velocity and density noise in the bulk flow. It is also important that the test posed here, and those commonly used in other works, are actually posed in a periodic domain with a compressible flow. To obtain a basis for comparison we use a numerical reference solution to the problem as posed and establish the uncertainty on this reference solution in a rigorous manner in Section 5.

To extract the amplitude of the y-velocity mode of the instability a more involved calculation is required. We wish to define the measurement in a manner that can be made consistent across different types of discretizations. A simple Fourier transform defined on a grid would be entirely appropriate for point-based finite-difference schemes or pseudospectral schemes, but is somewhat less well motivated for finite-volume schemes, and inappropriate for meshless or unstructured mesh schemes. To state the analysis in a manner that is straightforward to describe for all codes, which treats all results in the same manner we use a discrete convolution. For the case of a uniform grid this amplitude M is given by

$$\begin{eqnarray} s_i &=\left\lbrace \begin{array}{@{} ll}\! V_y\sin (4\pi x)\exp (-4\pi |y-0.25|) & \! \! {\rm if }\; y<0.5\\ \ms \! V_y\sin (4\pi x)\exp (-4\pi |(1-y)-0.25|) & \! \! {\rm if }\; y \ge 0.5\\ \end{array}\right. \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} c_i &=\left\lbrace \begin{array}{@{} ll}\! V_y\cos (4\pi x)\exp (-4\pi |y-0.25|) &\! \! {\rm if }\; y<0.5\\ \ms \! V_y\cos (4\pi x)\exp (-4\pi |(1-y)-0.25|) &\! \! {\rm if }\; y\ge 0.5\\ \end{array}\right. \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} d_i &=\left\lbrace \begin{array}{@{} ll}\exp (-4\pi |y-0.25|) & {\rm if }\; y<0.5\\ \ms \exp (-4\pi |(1-y)-0.25|) & {\rm if }\; y\ge 0.5\\ \end{array}\right. \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} M &= 2\sqrt{\left(\frac{\sum _{i=1}^N{s_i}}{\sum _{i=1}^N{d_i}}\right)^2+\left(\frac{\sum _{i=1}^N{c_i}}{\sum _{i=1}^N{d_i}}\right)^2}, \end{eqnarray} \tag{ 9 }$$

where i ranges over all grid points or cell centers (N total grid points or cell centers) and the positions (x, y) are grid points or cell centers. This expression can be used in two or three dimensions. The mode amplitude M is calculated at each time snapshot.

For SPH simulations, each particle needs a different weighting in the sums as the particle density varies. In the case of variable smoothing length SPH where the smoothing length is set to encompass a fixed number of neighbors, we can use the smoothing length h_i for particle i to do this weighting. In the following formulae, p is the number of dimensions the simulation is run in

$$\begin{eqnarray} s_i &=\left\lbrace \begin{array}{@{} ll}\! V_y h_i^p \sin (4\pi x)\exp (-4\pi |y-0.25|) &\!\! {\rm if }\; y<0.5 \nonumber\\ \ms \! V_y h_i^p \sin (4\pi x)\exp (-4\pi |(1-y)-0.25|) &\!\! {\rm if }\;y\ge 0.5 \end{array}\right.\\ \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} c_i &=\left\lbrace \begin{array}{@{} ll}\! V_y h_i^p \cos (4\pi x)\exp (-4\pi |y-0.25|) &\!\! {\rm if }\; y<0.5 \nonumber\\ \ms \! V_y h_i^p \cos (4\pi x)\exp (-4\pi |(1-y)-0.25|) &\!\! {\rm if }\;y\ge 0.5 \end{array}\right. \\ \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} d_i &=\left\lbrace \begin{array}{@{} ll}\! h_i^p \exp (-4\pi |y-0.25|) & {\rm if }\; y<0.5\\ \ms \! h_i^p \exp (-4\pi |(1-y)-0.25|) & {\rm if }\;y\ge 0.5\\ \end{array}\right. \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} M &= 2\sqrt{\left(\frac{\sum _{i=1}^N{s_i}}{\sum _{i=1}^N{d_i}}\right)^2+\left(\frac{\sum _{i=1}^N{c_i}}{\sum _{i=1}^N{d_i}}\right)^2}. \end{eqnarray} \tag{ 13 }$$

The quantities s_i and c_i are defined for each particle i = 1, ..., N, from the position (x, y) and the y-velocity V_y of that particle. An advantage of the definition used here is that we can directly analyze the SPH particle values as simulated without introducing an additional interpolation to a fixed grid. This feature follows over to a unstructured mesh or meshless code.

For an unstructured mesh code, or a meshless code that defines quadrature volumes for the points the appropriate general form would be

$$\begin{eqnarray} s_i &=\left\lbrace \begin{array}{@{} ll}\! V_y w_i \sin (4\pi x)\exp (-4\pi |y-0.25|) & {\rm if }\; y<0.5 \nonumber\\ \ms \! V_y w_i \sin (4\pi x)\exp (-4\pi |(1-y)-0.25|) & {\rm if }\;y\ge 0.5 \end{array}\right. \\ \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} c_i &=\left\lbrace \begin{array}{@{} ll}\! V_y w_i \cos (4\pi x)\exp (-4\pi |y-0.25|) & {\rm if }\; y<0.5 \nonumber\\ \ms \! V_y w_i \cos (4\pi x)\exp (-4\pi |(1-y)-0.25|) & {\rm if }\;y\ge 0.5 \end{array}\right. \\ \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} d_i &=\left\lbrace \begin{array}{@{} ll} \! w_i \exp (-4\pi |y-0.25|) & {\rm if }\; y<0.5\\ \ms \! w_i \exp (-4\pi |(1-y)-0.25|) & {\rm if }\; y\ge 0.5\\ \end{array}\right. \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} M &= 2\sqrt{\left(\frac{\sum _{i=1}^N{s_i}}{\sum _{i=1}^N{d_i}}\right)^2+\left(\frac{\sum _{i=1}^N{c_i}}{\sum _{i=1}^N{d_i}}\right)^2}, \end{eqnarray} \tag{ 17 }$$

where w_i is the area or volume of cell or the quadrature volume for point i and the positions are the cell centers or point positions.

For an infinite domain with incompressible flow, the growth of the velocity mode M should be ∝exp (4.384 × t) (Wang et al. 2010; Equation (18)). The growth rate for a KHI with these two conditions has been used before as a comparison for results obtained in periodic domains with a compressible flow, but the two problems are formally different, and depending on the parameters the growth rates may vary. Again, to circumvent this difficulty, we compare results to the numerical solution of the test problem specified, and establish the uncertainty on this reference solution in a rigorous manner in Section 5.

In this paper, we compare the results from six codes to the reference solution, which itself is produced with the Pencil code.

The Pencil code⁶ is a fixed Eulerian mesh, non-conservative, finite-difference, MHD code that uses sixth-order centered spatial derivatives and a third-order Runge–Kutta time-stepping scheme, being primarily designed for weakly compressible turbulent hydromagnetic flows. For the problem in question, in order to keep the Reynolds number low at the grid scale while keeping the integral and intermediate scales nearly inviscid, explicit sixth-order hyperdiffusion and hyperviscosity are added to the mass and momentum equations as specified in Lyra et al. (2008).⁷

The other codes, Enzo, Athena, NDSPMHD, and Phurbas are introduced below.

Enzo is a three-dimensional, Eulerian adaptive mesh refinement hybrid (hydrodynamics + N-body) grid-based code (Bryan & Norman 1995; O'Shea et al. 2004).⁸ For this problem, the Euler equations are solved using a third-order piecewise parabolic method (PPM) with the two-shock approximate Riemann solver. Time stepping is constrained by a Courant condition for the gas with a Courant factor C = 0.4. The run-time PPM diffusion, flattening, and steepening parameters were set to zero. Enzo version 1.5 was used.

Athena is a three-dimensional Eulerian grid code that (among other algorithms) implements a higher order Godunov method for hydrodynamics (Stone et al. 2008). Specifically, we have used the third-order cell reconstructions with the HLLC approximate Riemann solver and the unsplit corner–transport–upwind (second-order time integration algorithm. Otherwise the options used were as specified in the two-dimensional test problem supplied with the code, with a Courant number C = 0.8. We used Athena version 4.1 obtained from the project Web site.⁹

NDSPMHD is a one-, two-, and three-dimensional reference implementation of SPH and a platform for experimentation (Price 2012). We obtained NDSPMHD version 1.0.1 from the author's Web site.¹⁰ NDSPMHD was run on this problem in two dimensions, using both the cubic- and quintic-kernel options. The cubic kernel is the conventional choice for SPH, whereas the quintic kernel delivers higher accuracy at the cost of computational expense. Price (2012) describes the NDSPMHD implementation of SPH as converging as higher order kernels are used. That is, the result on the test problem shown here should converge with the combination of using more particles and using a higher order kernel. NDSPMHD also supports the artificial thermal conductivity described in Price (2008). The results of SPH simulations may depend strongly on the initial particle distribution used. To generate the initial condition, we relaxed a set of equal mass particles to an approximate equilibrium with an artificially imposed pressure field which produced the required density profile. The particles settle into a roughly hexagonal grid, although with dislocations required to produce the spatially varying density. The number of particles used at each resolution matched the number of cells or points used for the grid code (128², 256², 512²). Otherwise, the code was run with the default parameters used in Test 6 of the NDSPMHD examples package.

Phurbas is a meshless, adaptive, Lagrangian code for magnetohydrodynamics (Maron et al. 2011; McNally et al. 2011). Phurbas uses third-order least-square fits to derive spatial derivatives, and a second-order scheme for time integration. Stabilization is achieved through an artificial bulk viscosity. It is run here in three dimensions, using volumes with height 1/64, 1/128, and 1/256 in thickness in the z-direction. Phurbas does not use a grid, so instead we use spatially constant resolution and set the resolution parameter λ to the cell size used in the grid codes. To produce the initial particle distribution, we first used a tiling procedure as in McNally et al. (2011) and then further relaxed the distribution to one that would arise naturally in a shearing flow by running the problem to t = 1.5 and restarting the test with the initial condition defined on resulting particle distribution. Because the disordered particle distribution is inherently three-dimensional, the results at a given resolution cannot be strictly compared to the two-dimensional runs performed in other codes here.

To produce a solution to the full nonlinear, periodic, compressible case as run in this work, we performed an extensive convergence study with the Pencil Code. This convergence study allows us to establish not only a very high quality reference solution, but also a notion of the uncertainty in this reference solution. The importance of the unusual step of establishing the uncertainty of the reference result is that we can then assert with confidence that the differences seen between other lower quality results and this reference result are overwhelmingly due to errors in the lower quality solutions. In Section 6, the Pencil Code is shown to be well suited to the smooth, subsonic problem posed here. We use grids of 128 × 128, 256 × 256, 512 × 512, 1024 × 1024, 2048 × 2048, and 4096 × 4096 points, specified so that every second grid coordinate overlaps on successive refinements, and with the time-stepping scheme in the Pencil Code modified to provide outputs at exact Δt = 0.02 time unit intervals. This set of outputs enables a resolution study at each output time for the convergence of the mode amplitude. Establishing the empirical rate of convergence of the mode amplitude M allows a Richardson extrapolation-based estimate of the uncertainty in the most resolved measurement. Hence, we are able to make comparisons of the results from other codes to the highest resolution Pencil Code result while knowing in a rigorous manner that the errors in this reference result are negligible.

First, we can calculate the empirical rate of convergence p of the mode amplitudes defined by Equations (6)–(9) for a set of three results with a refinement ratio of 2 between each resolution as

$$\begin{equation} p = \ln \left(\frac{f_3-f_2}{f_2-f_1}\right)\bigg/ \ln (2), \end{equation} \tag{ 18 }$$

where f₃ is the value of the mode amplitude on the coarsest grid, and f₂ and f₁ are the values on the medium and finest grid, respectively (Roache 1998; Equation (5.10.6.1)). This rate of convergence tells us how fast the series of values from each resolution is converging toward the correct result. Once we have identified the convergence rate of the series of results, we can apply a generalized form of Richardson extrapolation to estimate the converged result, and hence derive an indication of the uncertainty in our highest resolution result. This indication of uncertainty is the Grid Convergence Index (GCI; Roache 1998), a uniform method of reporting the uncertainty on such a convergence study given as

$$\begin{equation} {\rm GCI} = F_s\frac{(f_2-f_1)/f_1}{1-r^p} \end{equation} \tag{ 19 }$$

from Roache (1998; Equation (5.6.1)). The value of the safety factor F_s we use is 1.25. This value is that suggested by Roache (1998; Section 5.9) as being appropriate when the rate of convergence is explicitly determined with a convergence study as in this work.

A density plot at time t = 1.5 from the highest resolution (4096²) calculation is shown in Figure 2. The results of evaluating Equation (18) for each set of three resolutions is shown in Figure 3. This figure shows that the convergence rate settles at approximately 1.75 for most of the time interval t = 0–1.5 when the highest resolution results are considered. Using the observed rate of convergence at each time, we can assign the uncertainty on the result with Equation (19), which is shown in Figure 4. The high-resolution results used are necessary to establish a well-behaved convergence in Figure 3, which means that the uncertainty is only well known when the uncertainty itself is very small. To demonstrate more explicitly the convergence behavior and the magnitude of the changes between successive resolutions, we have plotted the differences in the y-velocity values between successive resolutions in one quadrant of the domain in Figure 5. The greatest changes between successive resolutions are localized to the density change interface and show not suggestion of the presence of secondary instabilities. A similar plot for the density is shown in Figure 6, again showing no suggestion of secondary instabilities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulations are identified by a two letter prefix as outlined in Table 1 and the resolution (128², 256², 512²). Results for the y-velocity unstable mode amplitude and the growth of that quantity for all codes are plotted in Figure 7. Figure 8 gives the results for all codes for the minimal y-direction kinetic energy. The following two subsections discuss these two measured quantities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. y-velocity Unstable Mode Amplitude

6.2. Maximum y-direction Specific Kinetic Energy

6.3. Density at t = 1.5

We show gray-scale slices of the density field at t = 1.5 in Figure 9. All the images have the same limits on the gray scale between density of 0.9883 and 2.0320, the density extremes in the highest resolution result in the Pencil Code convergence study. The results for Pencil, Enzo, and Athena are largely similar, as at high resolution these codes agree well with each other and with the reference result. Though the result with Phurbas strongly resembles the reference result, it clearly shows more diffusion. The SPH results from NDSPMHD (only Ne and No shown) reflect the slow growth of the unstable y-velocity mode already discussed. The simulation No result using quintic-kernel SPH shows less diffusion than simulation Ne using cubic-kernel SPH. Especially in simulation Ne, secondary features of a filamentary appearance can be seen along the interface, and these are less apparent in the simulation No result. The quintic-kernel result (No) overall shows better agreement with the reference than the cubic-kernel result (Ne).

Zoom In Zoom Out Reset image size

Overall, the grid-based codes Pencil Code, Athena, and Enzo had very similar performances. For these codes, the test problem in this work (run to t = 1.5) confirms their correctness. This shows that the test as outlined here can be used to discriminate among numerical schemes. In this test, we demonstrated that Phurbas and NDSPMHD, while both using meshless Lagrangian schemes, give significantly different convergence behaviors. Though Phurbas was run in three dimensions, and NDSPMHD in two, the strikingly different qualitative behavior bears some explanation. A primary observation is that Phurbas differs from NDSPMHD in that Phurbas uses a third-order accurate and consistent spatial discretization, while NDSPMHD uses an SPH discretization which has zeroth-order inconsistency. This issue is sufficiently complex that it is discussed in a separate section (Section 8). We also note that no code developed obvious signs of secondary instabilities in the solution by time t = 1.5, in agreement with the findings of the convergence study performed on the reference result. How, and when, secondary instabilities may arise in a KHI test such as this is discussed in Section 9.

The results for the maximal y-direction kinetic energy in Section 6.2 show significant noise appearing in the velocity in the SPH simulations Ne, Nc, and No. This section is devoted to exploring the source and behavior of this noise. It has been argued that maintaining particle order is vital to achieving good results with SPH (Price 2011). Particle ordering in SPH can be expressed as a condition that the Lagrangian of the system of particles is minimized (Price 2011; Section 2.5). To seek this minimum, the particles must have some re-meshing motion in addition to the pure fluid motions (Price 2012; Section 5.2). These re-meshing motions mean that in SPH one always has some motions that are not physical, but purely related to the SPH particles attempting to relax to an ordered state (Price 2012; Section 5.2). These re-meshing motions are also shown in the post-shock state in Price (2012; Figure 10). The re-meshing motions are provided by the linear errors in the SPH pressure forces, which are in turn a result of the zeroth-order inconsistency of SPH interpolation. That is, the zeroth-order inconsistency in SPH interpolation provides a linear error in the pressure force which causes two particles approaching each other to repel. Re-meshing motions created by this repulsion are in turn damped by the artificial viscosity to encourage the particle distribution to relax. In this way, both the zeroth-order inconsistency in the pressure estimate is vital to maintaining particle order and the artificial viscosity cannot simply be disabled. Further, though more advanced artificial viscosities can be designed, the identification of the particle velocity with the fluid velocity means that the need for motions preserving particle order will necessarily corrupt, to some degree, the fluid velocity itself. The root cause that creates this situation is the zeroth-order inconsistency in SPH interpolation, so the parameter we vary is the one that varies this error, the choice of SPH smoothing kernel. As the change from the cubic kernel to the quintic kernel decreases the size of the zeroth order inconsistency, the re-meshing pressure forces are smaller, and the resulting velocities are smaller. Consequently, the level of y-direction kinetic energy noise seen in the simulation No is smaller than that seen in the simulation Ne.

Recently, Cha et al. (2010) and Read et al. (2010) have connected zeroth-order inconsistency in SPH to poor results for related KHI test problems, and Bauer & Springel (2012) has demonstrated the connection in the context of low Mach number turbulence. In this work, we have demonstrated this connection by first showing that on a smooth test problem the artificial conduction of Price (2008) does not significantly affect the results. Then, we have demonstrated that when the quintic kernel is used, with smaller inconsistency errors than the cubic kernel, the test results improve significantly. The size of the inconsistency errors are reflected in the maximal y-direction specific kinetic energy statistic, as pressure gradient errors dive spurious particle motion. Finally, Phurbas, while being meshless and Lagrangian-like SPH, uses a third-order consistent interpolation. The consequence of this for the test in this paper is that it performs much better than SPH, as the pressure forces are accurate enough to keep the velocity and density noise much smaller than in SPH. Hence, Phurbas has a qualitatively different behavior on this test, and convergence does not depend on varying the number of neighboring particles used in the interpolation.

To demonstrate that the velocity noise behavior seen in this work is general, we have run a series of additional tests. The first test is a version of the KHI setup in Section 2 with a uniform density of 1.0. For SPH, this is a particularly simple choice as a uniform hexagonal close packed grid of particles is the unique relaxed distribution in two dimensions. Hence, initially the setup does not suffer from any velocity noise. Figure 10 shows that regardless of this initially relaxed distribution, the maximal y-direction kinetic energy still reflects the growth of the velocity noise. Again, as in the previous tests the noise grows sooner at higher resolutions.

Zoom In Zoom Out Reset image size

We note that the velocity noise in the highest resolution quintic-kernel case of our isodensity test appears to be triggered by the growth of the primary KHI instability. To simplify the setup further, we remove the y-direction velocity perturbation from the initial condition, yielding a smooth unperturbed shear flow. For the maximum y-direction specific kinetic energy measurement, the trivial analytic solution for this problem is a value of zero for all times. Figure 11 (upper panel) shows that in this setup, run with the cubic kernel, the maximal y-direction kinetic energy grows to the same level as before in the isodensity KHI test, although it takes longer. This growth happens at earlier times for higher resolutions. Figure 11 (lower panel) displays the same behavior for the quintic kernel, although the timescales involved are longer. These simulations are stopped abruptly when they succumb to particle pairing instability (Price 2012; Section 5.4) and two particles approach within 10⁻⁸ length units. Some recent proposed modifications to SPH reduce or eliminate this instability for large kernels (Read et al. 2010; Read & Hayfield 2012; Dehnen & Aly 2012). With each kernel choice, the time interval until the velocity noise jumps is shorter as the number of particles is increased, but for a given number of particles the time interval until the velocity noise jumps is longer if the quintic kernel is used. That is, the results converge toward analytic solution as the zeroth-order inconsistency in the SPH interpolation is reduced.

Zoom In Zoom Out Reset image size

We have shown that given a convergence test stated in a well-posed manner, all the methods tested appear to converge toward the correct result for the growth of the primary instability. Recent discussion of Kelvin–Helmholtz tests has broadened to include secondary instabilities. Springel (2011) shows secondary instabilities developing from a similar initial condition. The reference solution we compute shows no indication of these structures. Springel (2011) suggested that their moving-mesh code is able to resolve secondary KHI billows which cannot be resolved in their fixed-mesh code because it was too diffusive. Though the solution in a fixed-mesh and moving-mesh code should not be expected to be equivalent at any finite resolution, that a given code does not develop secondary KHIs is not simply a function of the diffusivity of the code, it is also a dependent on the seeding of such instabilities.

In general, we can categorize the possibilities for why our reference case and the tests done at lower resolutions in this work do not show the development of any secondary instabilities in three cases.

• 1.  
  The secondary billows should grow physically, either due to the nature of the initial perturbation fed into the problem or the interaction between some combination of the initial perturbation and modes of the instability directly seeded by the initial perturbation. In this case, the secondary billows should eventually show up at some resolution in any convergent code, but should arise at a particular location and time.
• 2.  
  The secondary billows grow due to the balance of numerical perturbations and numerical diffusion. In this case the billows seen at some time should disappear at some resolution in any convergent scheme as the significant power in the numerical perturbations should eventually move to spatial scales too small to seed the secondary instability efficiently.
• 3.  
  The slight differences between the setup of Springel (2011) and our setup mean the difference between seeing physical growth of secondary billows and failing to produce them.

In the first case, developing the secondary instabilities is merely a matter of using a sufficiently large resolution. At lower resolutions, a resolution study should still suggest a significant uncertainty or change between simulations of different resolution as the secondary billows are damped less and less. In the second case, the resolution necessary to make the secondary billows disappear may be quite large, as a numerical mechanism introducing noise at a small scale may still have significant power at larger wavelengths depending on the spatial correlation of the mechanism introducing the noise. However, the resolution study performed with the Pencil Code to much higher resolution shows no indication that these modes grow. The third case cannot be ruled out explicitly, as Springel (2011) does not specify the exact details of the setup used. However, we can show that for our problem that secondary instabilities that do develop are of a purely numerical origin. This strongly suggests that the secondary billows seen in Springel (2011) are a numerical artifact, so the observation that a fixed grid codes does not develop these on the same problem does not imply that the fixed grid code is too diffusive to support the modes.

To demonstrate how numerical effects can seed secondary KHI we have performed a test with Athena. We ran the KHI test at resolutions of 1024², 2048², and 4096² until t = 4. In Figure 12, the density in a region centered on a single primary KHI billow is shown at time t = 3. Secondary KHI billows can be seen growing in the 1024² case, and this pattern is successively suppressed at the higher resolutions, suggesting it is an artifact of the finite resolution and is converging away. However, in the 2048² resolution, a different set of secondary instabilities can be seen growing at much shorter wavelengths in the central winding of the primary billow. As the resolution is increased, the numerical seeding of the secondary instabilities changes, and the secondary modes that are excited change. Figure 13 shows the same region at time t = 3.2. By this point, the secondary instabilities in the 4096² resolution simulation have become apparent. Surprisingly, a new set of secondary billows has appeared on the outer winding of the primary billow. We cannot reach well-justified conclusions about a particular mode of the secondary instability from this study as we cannot reproduce the same instability at two different resolutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Though the growth of secondary instabilities is likely a physical reality at the Reynolds numbers involved in astrophysical problems, relying on numerical effects to seed them will not result in a true physical model of the phenomena as the seeding, and hence growth of these instabilities, will be inherently dependent on parameters such as resolution. Models relying on numerically seeded instabilities, even if the presence of the instability is physical, make it difficult to separate numerical effects from physical behavior, which in turn makes it difficult to come to strong conclusions about the effect of the instability. Conclusions about the instability must consist of a measurement and some manner of characterization of the error in that measurement. In order to characterize the error in the model of the instability, a convergence study must be performed. To perform this convergence study, a fixed seeding of the instability must be possible across all resolutions. If the seeding is a numerical effect caused by the finite resolution, it will not be fixed between two resolutions. Hence, the tests in this work do not show that any code used cannot, at any of the resolutions tested, resolve secondary KHIs as these have not been seeded in a controlled way or even in an avoidable manner.

From the results of the convergence study, in Section 5 we propose that if a code develops secondary Kelvin–Helmholtz billows in this test by t = 1.5, it is due to the growth of numerical perturbations. The less rigorous study performed with Athena suggests that the same conclusion should hold to at least t = 2.5. In the limit of infinite resolution, any convergent code should reproduce the correct result. However, if at finite resolution a code shows a tendency to produce secondary instabilities, then the scheme can be improved by adding a diffusive operator to damp the noise leading to the instability. Particularly with respect to moving-mesh tessellation codes, Kelvin–Helmholtz tests are not the only place where behavior suggests that some additional numerical diffusion should be used to damp grid scale noise.

Development of secondary KHI after t = 1.5 is likely due to the presence of spurious noise in the solution. For example, in the preparation of this work, we discovered that the evolution of the test problem here differed greatly at high resolution (4096²) between Enzo versions 1.5 and 2.0. In Enzo 2.0, a bug existed that caused slightly incorrect pressure reconstructions. This caused small sound waves to launch from the interface, and propagate though the periodic domain interacting with themselves and forming small short-wavelength perturbations. The discovery of this bug was fortuitous, however, because it demonstrates again how artificial, numerical perturbations can give rise to secondary Kelvin–Helmholtz in this test problem if they are able to overwhelm the dissipation of the scheme.

The underlying cause of the tendency of many schemes to develop secondary KHI in this test problem is that the shear interface becomes increasingly steep as it stretches in the primary KHI billow. Eventually, the width of the interface approaches the grid scale and it become susceptible to numerically seeded secondary instabilities. This behavior is also commonly seen in the initial evolution of grid-aligned sharp transition versions of the KHI test. Two examples are Junk et al. (2010; Figure 13) and Springel (2011; Figure 8, upper right panel). We suggest then that even fixed-grid finite-volume Godunov schemes may be improved in pathological cases of unresolved shear interfaces by the addition of a diffusive flux. This flux should be chosen to spread the interface over enough grid cells to suppress the numerically seeded instabilities.

Another lesson to be derived here is a cautionary one. Not all new instabilities seen as resolution is increased when solving the discretized Euler equations are physically real. New numerical instabilities can reveal themselves as resolution is increased, as the flow can enter into new regimes where it is more sensitive to the inevitable numerical noise in a method. In the high-resolution set of Athena simulations, this can be seen in the magnitude of the density gradients. In Figure 14, the density gradient in a slice through a primary billow at time t = 3.0 in the three Athena simulations used in this section is plotted, calculated with a four point second-order finite-difference stencil. As the resolution is increased, the maximum gradient achieved increases. Mathematically, when solving the Euler equations, this behavior arises because the modified equations that are actually solved by the method change as resolution is increased—the diffusive effects become smaller. One route around this difficulty can be to solve the Navier–Stokes or Boltzmann equations instead with a fixed viscosity or particle mean free path. Since these equations have a physical scale where diffusion dominates dynamics, the reliable elimination of numerically generated instabilities for arbitrarily long run times can be obtained by fixing the physical diffusive scale and reducing the grid scale far below the diffusive scale.

Zoom In Zoom Out Reset image size

The same point illustrates how the transition to turbulence and mixing must be studied when the Euler equations are used. To produce the secondary instabilities that break up the flow, the nonlinear interaction or modes or seeding of secondary instabilities must be done on a controlled manner. This job cannot be left to the numerical noise, or the time and manner in that the flow breaks up will be a reflection of numerical issues and not of physical reality.

Finally, we suggest that it is possible to produce a controlled test of the growth of secondary billows from definite perturbations, similar to the study performed by Fontane & Joly (2008) in an incompressible flow. Such a setup could be useful in determining the appropriate and minimal diffusion to add to a scheme to suppress the numerical seeding of secondary instabilities in given conditions.

We have constructed a reference solution with a well-characterized uncertainty, along with defining a general manner in that the test can be analyzed. This methodology was applied to example codes from the major families of numerical techniques used in astrophysics. All codes tested showed convergence toward the reference result when the resolution was increased in the appropriate manner. For SPH, the use of an artificial thermal conductivity does not significantly effect the results, but using a higher-order kernel (and hence a larger number of neighbors) does improve the results. We conclude then that the fundamental reason for poor performance of SPH in KHI is the zeroth-order inconsistency of SPH interpolation. Visually, to time t = 1.5 in the test problem there are no secondary instabilities that arise in the reference solution. By examining the relative behavior of different types of code, we argue that the presence of secondary instability on this test is caused by having a numerical diffusion that is very low compared to the grid noise in the method. Hence, we propose that it is advantageous in some methods, particularly moving-mesh tessellation methods, but also in fixed-grid Godunov schemes, to include an extra diffusion operator to smooth the solution such that grid noise does not drive small-scale instabilities.

We are indebted to the authors of NDSPMHD, Athena, Enzo, and the Pencil Code for making the codes freely available to us. The SPH visualizations were created with SPLASH (Price 2007). We thank the anonymous referee for constructive comments that significantly improved the organization of the paper. We thank Daniel Price for feedback on a draft of this manuscript, and suggesting the use of fully relaxed SPH initial conditions and the isodensity special case. We thank Mordecai-Mark Mac Low for his support and useful discussions. We thank Paul Duffel for useful discussions about moving-mesh schemes. J-C.P. thanks Orsola De Marco and Falk Herwig for their support. This work has been supported by National Science Foundation (NSF) grants AST-0835734 and AST-0607111. W.L. gratefully acknowledges partial financial support from the NSF under grant No. AST-1009802. W.L. completed co-writing this work at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. OCI-1053575. This work used computer time provided partially by Westgrid and Compute Canada.