{'fraction_non_alphanumeric': 0.034941329856584095, 'fraction_numerical': 0.019991308126901346, 'mean_word_length': 4.468631178707224, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical simulations of the Kelvin-Helmholtz instability with the Gadget-2 SPH code

The method of Smoothed Particle Hydrodynamics (SPH) has been widely studied and implemented for a large variety of problems, ranging from astrophysics to fluid dynamics and elasticity problems in solids. However, the method is known to have several deficiencies and discrepancies in comparison with traditional mesh-based codes. In particular, there has been a discussion about its ability to reproduce the Kelvin-Helmholtz Instability in shearing flows. Several authors reported that they were able to reproduce correctly the instability by introducing some improvements to the algorithm. In this contribution, we compare the results of Read et al. (2010) implementation of the SPH algorithm with the original Gadget-2 N-body/SPH code.


Introduction
The Kelvin-Helmholtz Instability (KHI) is the instability that appears at the interface between two shearing fluid flows of different densities. Many experimental and numerical results have been published where such instability is reproduced. In particular, in astrophysics such instabilities may be responsible for many phenomena observed in regions of high gas density contrast (Murray et al., 1993). The necessity for subsonic velocities for the KHI survival was first explored in the context of astrophysical flows by Vietri, Ferrara & Miniati (1997). As discussed by these authors, the growth rate of the KHI is smaller than the gas speed of sound. The correct modeling of KHI is essential since the vorticity and shear flows appear in diverse hydrodynamic processes, such as for example, the onset of turbulence.
The Smoothed Particle Hydrodynamics method (SPH) is a Lagrangian meshless particle method used for simulation of fluid transport. Since its original formulation the SPH method has been constantly improved. In recent years a search for possible differences between grid-based and particle-based methods has been widely discussed. For instance, Agertz et al. (2007) compare the results for several hydrodynamic tests obtained with SPH and grid-based methods. They found a striking difference between the results. In particular, the inability to reproduce the vorticity rolls in the shearing flow was claimed to be a deficiency of the SPH method. Price (2008) showed that introducing an artificial thermal conductivity term into the standard SPH, in order to smooth the discontinuities in the thermal energy, allowed for similar results to those obtained by Agertz et al. (2007) using grid-based simulations.
On the other hand, several authors suggested that the artificial conductivity is not the main factor responsible for suppression of the instability. It is well known that the choice of the smoothing kernel also affects the results. For example, the standard cubic spline kernel tends to suffer from the clumping instability (known also as pairing instability) for large number of neighbors, introducing errors and lowering the resolution. Read, Hayfield & Agertz (2010) showed that kernels that produce a constant force term in the center prevent the clumping. Several alternative kernel shapes were proposed to remedy this problem (see for example, Dehnen & Aly (2012)). Hubber, Falle, & Goodwin (2013) performed a comparison of KHI simulations obtained with SPH and with the AMR Eulerian code Pencil. They concluded that convergence between SPH and grid codes may be obtained if higher order kernels (i.e., quintic) that support larger numbers of neighbors are used. Cha, Inutsuka & Nayakshin (2010) showed that inaccurate density gradients that are obtained with the standard SPH formulation are responsible for the suppression of instabilities and can be alleviated by using a Godunov formulation of SPH. These works concluded that the standard SPH formulation is unable to reproduce the KHI and that some improvements must be implemented.
Some of the causes that suppress the KHI are the following: when pressure discontinuities appear as a result of the lack of entropy mixing on the kernel scale, when large errors are introduced in the momentum equation due to a finite number of neighbor particles, due to the pairing instability, or due to contact discontinuities.
In this work, we compare the results for the KHI as obtained using the standard Gadget-2 SPH with the formulation proposed by Read, Hayfield & Agertz (2010).

Initial conditions
As the initial conditions for the first part of the KHI tests (A), we use the example already included in the OSPH code (Read, Hayfield & Agertz, 2010). This is done with the aim of providing direct comparison of our simulations with previously published results. In the second part (B), we use the so-called "well-posed" initial conditions described in McNally, Lyra & Passy (2012) (see also Robertson et al. (2010)). The latter is regularized by smoothing the density and velocity on the interfaces in order to prevent abrupt pressure jumps. Such situations are frequently found for example in simulations of real astrophysical systems. In addition, the two sets of initial conditions differ by the type of density sampling, which is done by the spatial distribution of the particles (set A) or by varying the particle masses (set B). Set C differs from set B only in the magnitude of the densities and velocities.
For brevity, we describe only the initial conditions for the simulation sets B and C, which consist of a 3D slab of size 1.0 × 1.0 × 0.0325, with the central part having a density ρ 2 = 2ρ 1 and moving with the velocity v 2 = −u, while the upper and lower sections have a density ρ 1 and move with the velocity v 1 = u along the x-axis. A small velocity perturbation is added, v y = δ v y sin(2πx/λ ), with λ = 0.5. The gas has a constant pressure P = 2.5 everywhere and satisfies an ideal gas equation of state P = (γ − 1)Uρ, with γ = 5/3. The system of units is such that length, time and mass are equal to unity. For the simulation set B, we use ρ 1 = 1.0, u = 0.5 and δ v y = 0.01, while for the set C we define ρ 1 = 32.0, u = 0.1, δ v y = 0.01 and δ v y = 0.002. The system is sampled with an equally spaced cubic grid of 256 × 256 × 8 particles. The parameters are chosen with the purpose of staying in the subsonic regime, which guarantees the KHI formation.
The characteristic onset time of the KHI in the linear regime is given by: For the above initial parameters the characteristic times are τ A ≈ 3.4, τ B ≈ 1.06 and τ C ≈ 5.3, and the total simulation times were 4, 2.1, and 8 respectively. Note that different units are used for set A. We use the Gadget-2 code (Springel, 2005) and its modification -OSPH -the Optimized SPH introduced by Read, Hayfield & Agertz (2010). In both codes the flags -DPERIODIC, -DNOGRAVITY and -DLONG were activated, and in OSPH the recommended flags were also switched on (-DOSPHM, -DOSPHRT, and -DOSPHHOCT). The artificial viscosity was set to α = 0.8 in all simulations.

Results
We performed a number of tests using different combinations of the initial conditions and code parameters. In Table 1 we summarize some of them. The first column is the set of models coded by the first letter of the name, which are A, B or C, the second column gives the number of neighbors, and the third column shows the code employed for the simulation. Figure 1 (set A), Fig. 2 (set B) and Fig. 3 (set C) summarize the results of using different initial conditions and different numbers of neighbor particles. The runs with a higher number of neighbors with the standard Gadget-2 were not completed because a maximum number of tree-nodes was reached due to the pairing instability.
From Fig. 1, we observe that model AOSPH442 is the only run that produces well distinguished rolls, while using either a smaller number of neighbors or the Gadget-2 code leads only to mild perturbations. This fact confirms the results of Read, Hayfield & Agertz (2010). A notable feature is that in model AOSPH96 with a smaller number of neighbors the perturbations are completely damped. On the other hand, run AGADG64, which uses a nearly optimum number of neighbors for the 3D cubic spline kernel, shows pronounced undulations, although without any rolls. Using their new SPHS code, which implements a dissipation switch, Read & Hayfield (2012) obtained KHI in both single mass and multimass particle models. Their initial conditions for single mass particles are identical to those used in set A. Note that in the set A a Mach number for the low density layer is M ≈ 0.11.
The mass varying initial conditions in the simulations of set B lead to pronounced perturbations with well developed rolls being produced only in the runs with a maximum number of neighbors (Fig. 2). Small undulations observed at time t = τ B (left column) are evolved into a KHI for t = 2τ B (right column) in runs BGADG64 and BOSPH442. Surprisingly, standard SPH with multimass setup do develop KHI, con- trary to previous claims. Comparing a BGADG64 density projection at t = 2.1 to a reference solution obtained with the Pencil code of McNally, Lyra & Passy (2012) (their Fig. 2) we observe a very similar shape and amplitude of density rolls. However, the reference solution shows the density projection at t = 1.5, indicating that in BGADG64 the KHI develops slowly. In the case of BOSPH442, the density projection which matches approximately the reference solution corresponds to t = 1.7. In order to check the effect of a low number of neighbors we repeated the BGADG64 with 33 neighbors (not shown), and obtained a similar situation to BOSPH64 with no KHI. The Mach number of the low density layer is M 2 ≈ 0.34.
In order to explore the effect of different shear velocities and perturbation amplitudes, we performed some additional simulations with both codes reducing the initial velocities, and either reducing or keeping the same perturbation amplitude δ v y . These models are listed in Table 1 as set C, where the last letter in the name stands for large amplitude (L), δ v y = 0.01, and for small amplitude (S), δ v y = 0.002, respectively. In this case the gas is still subsonic with c 1 = 0.36 and c 2 = 0.25, giving the Mach number of low density layer M 2 ≈ 0.4. The Fig. 3 compares the results obtained with both codes. While these appear very similar, the amplitude of the KHI looks smaller and lacks secondary ripples for the Gadget-2 code.

Concluding remarks
The aim of this work was to test several implementations of the SPH code and determine the necessary conditions for successful simulations of the mixing process in shearing flows. We have compared visually the projected densities for standard and optimized SPH codes for non-linear structure formation. Runs with single mass standard SPH formulations showed a similar behavior to that described in previous studies, (c.f. Fig. 13 of Agertz et al. (2007) and Fig. 7 of Price (2008)). On the other hand, it seems that the multimass simulations do not depend on the SPH code implementation, but rather on the number of neighbors. The optimized SPH code is almost two times slower than the standard SPH for the same number of neighbors, -a disadvantage that was reported to be absent in the new SPHS code of the same authors (Read & Hayfield, 2012). We have shown that for mixing problems using the SPH formalism, it is essential to take much care in setting the initial conditions and the number of neighbors. Further investigation of the problem is necessary in order to obtain robust results with SPH.