Free Access
Issue
A&A
Volume 546, October 2012
Article Number A45
Number of page(s) 25
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201219715
Published online 02 October 2012

© ESO, 2012

1. Introduction

Smoothed particle hydrodynamics (SPH) is a Lagrangian, mesh-free particle method that is used to model fluid hydrodynamics in numerical simulations. The technique was originally developed in an astrophysical context (Lucy 1977; Gingold & Monaghan 1977), but since then it has been widely applied in many other areas (Monaghan 2005) of computational fluid dynamics.

The method has several properties that make its use particularly advantageous in astrophysical problems (Hernquist & Katz 1989; Rosswog 2009; Springel 2010a). Because of its Lagrangian nature, the development of high matter concentrations in collapse problems is followed naturally. Moreover, the method is free of both advection errors and Galilean invariant. Finally, the method naturally incorporates self-gravity and possesses very good conservation properties, its fluid equations being derived from variational principles.

The other computational fluid dynamical method commonly employed in numerical astrophysics is the standard Eulerian grid-based approach in which the fluid is evolved on a discretized mesh (Stone & Norman 1992; Ryu et al. 1993; Norman & Bryan 1999; Fryxell et al. 2000; Teyssier 2002). The spatial resolution of an Eulerian scheme based on a fixed Cartesian grid is often insufficient, however, to adequately resolve the high dynamic range frequently encountered in many astrophysical problems, such as galaxy formation. This has motivated the development of adaptative mesh refinement (AMR) methods, in which the spatial resolution of the grid is locally refined according to some selection criterion (Berger & Colella 1989; Kravtsov, Klypin & Khokhlov 1997; Norman 2005). Additionally, the order of the numerical scheme has been improved by adopting the parabolic piecewise method (PPM) of Colella & Woodward (1984), such as in the AMR Eulerian codes ENZO (Norman & Bryan 1999) and FLASH (Fryxell et al. 2000).

Application of these different types of hydrodynamical codes to the same test problem with identical initial conditions should in principle lead to similar results. However, there has been growing evidence over recent years that in a variety of hydrodynamical test cases there are significant differences between the results produced from the two types of methods (O’Shea et al. 2005; Agertz et al. 2007; Wadsley et al. 2008; Tasker et al. 2008; Mitchell et al. 2009; Read et al. 2010; Valcke et al. 2010; Junk et al. 2010).

For instance, Agertz et al. (2007) show that in the standard SPH formulations the growth of Kelvin-Helmholtz (KH) instabilities in shear flows is artificially suppressed when steep density gradients are present at the contact discontinuities. Moreover, the level of central entropy produced in binary-merger non-radiative simulations of galaxy clusters is significantly higher by a factor  ~2 in Eulerian simulations than in those made using SPH codes (Mitchell et al. 2009). The origin of these discrepancies has been recognized (Agertz et al. 2007; Read et al. 2010) as partly due to the intrinsic difficulty for SPH in properly modeling density gradients across fluid interfaces, which in turn implies that there is a surface tension effect that inhibits the growth of KH instabilities. A second problem is the Lagrangian nature of SPH codes, which prevents the mixing of fluid elements at the particle level and leads to entropy generation (Wadsley et al. 2008; Mitchell et al. 2009). In particular, Mitchell et al. (2009) argue that the main explanation of the different levels of central entropy found in cluster simulations is the different degree of entropy mixing present in the two codes. Unlike SPH, in Eulerian codes fluid mixing occurs by definition at the cell level and a certain degree of overmixing is certainly present in simulations made using mesh-based codes (Springel 2010b).

Given the advantages of the SPH codes highlighted previously, it appears worth pursuing any improvement in the SPH method capable of overcoming its present limitations. Along this line of investigation, many efforts have been made by a number of authors (Abell 2011; Price 2008; Wadsley et al. 2008; Valcke et al. 2010; Read et al. 2010; Heß & Springel 2010; Cha et al. 2010; Murante et al. 2011).

Price (2008) introduces an artificial heat conduction term in the SPH equations with the purpose of smoothing thermal energy at fluid interfaces. This artificial conductivity (AC) term in turn gives a smooth entropy transition at contact discontinuities with the effect of enforcing pressure continuity and removing the artificial surface-tension effect that inhibits the growth of KH instabilities at fluid interfaces. Similarly, Wadsley et al. (2008) suggest that in SPH the lack of mixing at the particle level can be alleviated by adding a heat diffusion term to the equations so as to mimic the effects of subgrid turbulence, thereby improving the amount of mixing.

Read et al. (2010) present an SPH implementation in which a modified density estimate is adopted (Ritchie & Thomas 1991), together with the use of a peaked kernel and a much larger number of neighbors. The authors show that the new scheme is capable of following the development of fluid instabilities in a more improved way than in standard SPH. Abell (2011) presents an alternative derivation of the SPH force equation that avoids the problem encountered by standard SPH in handling fluid instabilities, although the approach is inherently neither energy nor momentum conserving and is prone to large integration errors when the simulation resolution is low.

The method proposed by Inutsuka (2002) reformulates the SPH equations by introducing a kernel convolution so as to consistently calculate the density and hydrodynamic forces. The latter are determined using a Riemann solver (Godunov SPH). The method was revisited by Cha et al. (2010) and Murante et al. (2011), who show that the code correctly follows the development of fluid instabilities in a variety of hydrodynamic tests.

A deeper modification than those presented here so far was introduced by Heß & Springel (2010), who replaced the traditional SPH kernel approach with a new density estimate based on Voronoi tesselation. The authors show that the method is free of surface tension effects and therefore the growth rate of fluid instabilities is not as adversely affected as in standard SPH.

Finally, a radically new numerical scheme was developed by Springel (2010b) to retain the advantages of both SPH and mesh-based codes. In the new code, the hydrodynamic equations are solved on a moving unstructured mesh using a Godunov method with an exact Riemann solver. The mesh is defined by the Voronoi tesselation of a set of discrete points and allowed to move freely with the fluid. The method is therefore adaptative in nature and thus Galilean invariant but, at the same time, the accuracy with which shocks and contact discontinuities are described is that of an Eulerian code. Bauer & Springel (2012) argue that the standard formulation of SPH fails to accurately resolve the development of turbulence in the subsonic regime (but see Price 2012b, for a different viewpoint). The authors draw their conclusions by analyzing results from simulations of driven subsonic turbulence made using the new moving-mesh code, named AREPO, and a standard SPH code.

Similar conclusions were reached in a set of companion papers (Sijacki et al. 2012; Vogelsberger et al. 2011), in which the new code was used in galaxy formation studies to demonstrate its superiority over standard SPH. However, the code is characterized by considerable complexity making the use of the SPH scheme still appealing, and more generally, it is desirable that simulation results produced with a specific code should be reproduced with a completely independent numerical scheme when complex non-linear phenomena are involved. It appears worthwhile, therefore, to investigate, along the line of previous authors, the possibility of constructing a numerical scheme based on the traditional SPH formulation that is capable of correctly describing the development of fluid instabilities and at the same time incorporates the effects of self-gravity when present.

This is the aim of the present study, in which the SPH scheme is modified by incorporating an AC diffusion term into the equations as described by Price (2008). However, in Price (2008) the strength of the AC is governed by a signal velocity that is based on pressure discontinuities. For simulations where gravity is present, this approach is not applicable because hydrostatic equilibrium requires pressure gradients. An appropriate signal velocity for conductivity when gravity is present is then used to construct an AC-SPH code with the purpose of treating the growth of fluid instabilities self-consistently. The viability of the approach is tested using a suite of test problems in which results obtained using the new code are contrasted with the corresponding ones produced in previous work using different schemes. The code is very similar in form to that presented by Wadsley et al. (2008), but here the energy diffusion equation is implemented in a different manner.

The paper is organized as follows. Section 2 presents the hydrodynamical method and introduces the AC approach. In Sect. 3, we investigate the effectiveness of the method by presenting results from a suite of purely hydrodynamical test problems, namely the two-dimensional Kelvin-Helmholtz instability, the Sod shock tube, the point explosion or Sedov blast-wave test, and the blob test. In Sect. 4, we then discuss the results of hydrodynamic tests that include self-gravity. Specifically, we consider the cold gas sphere or Evrard collapse test, the Rayleigh-Taylor instability, and the cosmological integration of galaxy clusters. Finally, the main results are summarized in Sect. 5.

2. The hydrodynamic method

Here we present the basic features of the method and for reviews refer to Rosswog (2009), Springel (2010a), and Price (2012a).

2.1. Basic SPH equations

In the SPH method, the fluid is described by a set of N particles with mass mi, velocity vi, density ρi, and a thermodynamic variable such as the specific thermal energy ui or the entropic function Ai. The particle pressure is then defined as , where γ = 5/3 for a monoatomic gas. The density estimate ρ(r) at the particle position ri is given by (1)where rij ≡ ri − rj, W( | rij | ,hi) is the interpolating kernel that has compact support and is zero for  | rij |  ≥ ζhi (Price 2012a). The kernel is normalized to the condition WdDr = 1. The sum in Eq. (1) is over a finite number of particles and the smoothing length hi is a variable that is implicitly defined by the equation (2)where D is the number of spatial dimensions, fD = π, 4π/3 for D = 2, 3, and Nsph is the number of neighboring particles within a radius ζhi. This equation can be rewritten by defining hi in units of the mean interparticle separation (3)so that and . The smoothing length hi is determined by solving the non-linear Eq. (2), and it should be kept in mind that Nsph does not necessarily need to be an integer but can take arbitrary values if η is used as the fundamental parameter determining hi (Price 2012a). A kernel commonly employed is the M4 (cubic spline), which is zero for ζ ≥ 2. In 3D, typical choices of Nsph lie in the range Nsph ~ 33 − 64, which for the M4 kernel corresponds to η ≃ 1. − 1.25. The equation of motion for the SPH particles can be derived from the Lagrangian of the system (Springel 2010a; Price 2012a) and is given by (4)where the coefficients Ωi are defined as (5)These terms are present in the momentum Eq. (4) because the smoothing length hi itself is implicitly a function of the particle coordinates through Eq. (2).

2.2. Artificial viscosity

In SPH, the momentum equation must be generalized by adding an appropriate viscous force, which is aimed at correctly treating the effects of shocks. An artificial viscosity (AV) term is then introduced with the purpose of dissipating local velocities and preventing particle interpenetration at the shock locations. The new term is given by (6)where the term is the symmetrized kernel and Πij is the AV tensor. In the SPH entropy formulation (Springel & Hernquist 2002), it is the entropy function per particle Ai that is integrated and its time derivative is calculated as follows (7)where vij = vi − vj. For the AV tensor, we adopt here the form proposed by Monaghan (1997) in analogy with the Riemann solvers (8)where (9)is the signal velocity and μij = vij·rij/ | rij |  if vij·rij < 0 but zero otherwise, so that Πij is non-zero only for approaching particles. Here scalar quantities with the subscripts i and j denote arithmetic averages, ci is the sound speed of particle i, the parameter αi regulates the amount of AV, and fi is a controling factor that reduces the strength of AV in the presence of shear flows. The latter is given by Balsara (1995), as (10)where (·v)i and ( × v)i are the standard SPH estimates for divergence and curl (Monaghan 2005). For pure shear flows,  |  × v | i ≫  | ·v | i, so that the AV is strongly damped.

Using the signal velocity (9), the Courant condition on the timestep of particle i is given by (11)In the standard SPH formulation, the viscosity parameter αi controling the strength of the AV is given by αi = constant ≡ α0, with α0 = 1 being a common choice (Monaghan 2005). This scheme has the disadvantage that it generates viscous dissipation in regions of the flow that are not undergoing shocks. To reduce any spurious viscosity effects, Morris & Monaghan (1997) proposed making the viscous coefficient αi time-dependent so that it can evolve in time under certain local conditions. The authors proposed an equation of the form (12)where Si is a source term and τi regulates the decay of αi to the floor value αmin away from shocks. For the source term, the adopted expression is (13)which is constructed in such a way that it increases in the presence of shocks. The prefactor S0 is unity for γ = 5/3 and the damping factor fi is inserted to account for the presence of vorticity. The original form proposed by Morris & Monaghan (1997) was refined by the factor (αmax − αi) (Rosswog et al. 2000), which has the advantage of being more sensitive to shocks than the original formulation and of preventing αi from becoming higher than a prescribed value αmax. Cullen & Dehnen (2011) presented an improved version of the Morris & Monaghan (1997) AV scheme, employing as a shock indicator a switch based on the time derivative of the velocity divergence.

The decay parameter τi is of the form (14)where ld is a dimensionless parameter that controls the decay timescale. In a number of test simulations, Rosswog et al. (2000) found that appropriate values for the parameters αmax,αmin, and ld are 1.5, 0.05, and 0.2, respectively. In principle, the effects of numerical viscosity in regions away from shocks can be reduced by setting ld to higher values than ld = 0.2. In practice, the minimum time necessary to propagate through the resolution length hi sets the upper limit ld = 1. However, the time evolution of the viscosity parameter can be affected if very short damping timescales are imposed. Neglecting variations in the coefficients, the solution to Eq. (12) at times t > tin can be written as (15)where (16)and qi is a modified source term (17)From Eq. (15), it can be seen that αi(t) ≃ αmax in the strong shock regime Siτi ≫ 1 but this condition is not satisfied if Siτi ≲ 1. Therefore, for mild shocks this implies that, if very short decay timescales are imposed, the peak value of αi(t) at the shock front might be below the AV strength necessary to properly treat shocks.

To overcome this difficulty, a modification to Eq. (13) was adopted (Valdarnini 2011, hereafter V11) that, when ld = 1, compensates for the smaller values of Siτi with respect to the small ld regimes. This is equivalent to considering a higher value for αmax, so that in Eq. (13) αmax is substituted by αmax → ξαmax.

The correction factor ξ was calibrated using the shock-tube problem as reference and requiring a peak value of  ~0.6−0.7 for the viscosity parameter at the shock front, as in the ld = 0.2 case. The results indicate that a normalization of the form ξ = (ld/0.2)0.8 for ld ≥ 0.2 satisfies these constraints. This normalization has been validated in a number of other test problems showing that it is able to produce a peak value of the viscosity parameter at the shock location that is independent of the chosen value of the decay parameter ld.

The time-dependent AV formulation of SPH has been shown to be effective in reducing the damping of turbulence due to the effects of numerical viscosity in simulations of galaxy clusters (Dolag et al. 2005, V11). Moreover, it was used in Price (2012b) to argue that the conclusion of Bauer & Springel (2012) about the difficulty for SPH codes in properly modeling the development of subsonic turbulence is based on using only a SPH code in its standard AV formulation. In the following, unless otherwise stated, we describe simulations performed using a time-dependent AV, this will be fully specified by the set of parameters  { αmin,αmax,ld } .

2.3. Artificial conductivity

As already outlined in the Introduction, different formulations have been proposed for overcoming the problems encountered by standard SPH in the treatment of fluid discontinuities. Here we follow the approach suggested by Price (2008), who proposed adding a dissipative term to the thermal energy equation for smoothing the energy across contact discontinuities. The presence of these dissipative terms introduces a smoothing of entropy at fluid interfaces with the effect of removing pressure discontinuities.

The motivations behind this approach were discussed by Monaghan (1997), who showed how the momentum and energy equations in SPH must contain a dissipative term related to jumps in the variables across characteristics, in analogy with the corresponding Riemann solutions.

The artificial conductivity (AC) term for the dissipation of energy takes the form (18)where is the signal velocity, which does not need to be the same as that used in the momentum equation, eij ≡ rij/rij, and  is the AC parameter, which is of order unity1. Equation (18) represents the SPH analogue of a diffusion equation of the form (19)where (20)is the SPH expression for the Laplacian (Brookshaw 1985) and, in analogy with the analysis of Lodato & Price (2010) for defining an equivalent physical viscosity coefficient for the SPH numerical viscosity, we can define  as a numerical heat-diffusion coefficient given by (21)An important issue concerns the choice of the AC parameter , which must be constructed so that diffusion of thermal energy is kept under control and is damped away from discontinuities. This can be achieved by introducing a switch similar to that devised for AV. The dissipation parameter is then evolved according to (22)the meaning of the terms being similar to that of Eq. (12). In particular, for the decay timescale we here set and the floor value  to zero. For the source term , we use the expression (23)where fC is a dimensionless parameter of order unity and the parameter ε avoids divergences when ui → 0. This equation differs from the corresponding source term proposed by Price & Monaghan (2005) and Price (2008) by the factor , which has been inserted here in analogy with Eq. (13) for introducing a stronger response by the switch in the presence of discontinuities. The choice of values for the parameters fC and  depends on the problem under consideration, for example Price & Monaghan (2005) proposed fC = 0.1. Here a number of numerical experiments has shown that the best results in terms of the response of the AC parameter  to the presence of discontinuities are obtained by setting fC = 1 and , which we assume henceforth as reference values. Moreover, it is found that significant benefits in terms of the sharpness of the AC parameter profile across the discontinuity are obtained by inserting the term. In principle, the choice of the derivative term used to detect discontinuities is arbitrary, but in practice a second derivative (Price & Monaghan 2005) term ensures higher sensitivity to sharp discontinuities in thermal energy than a first derivative (Price 2005a).

An example of a signal velocity specifically designed to remove pressure gradients at contact discontinuities is that originally introduced by Price (2008) for pure hydrodynamical simulations (24)and further refined by Valcke et al. (2010). The ability of the new AC formulation of SPH to follow the development of KH instabilities using this expression for the signal velocity was verified in a number of tests (Price 2008; Valcke et al. 2010; Merlin et al. 2010). However, the disadvantage of using the signal speed (24) is that it cannot be applied when self-gravity is considered because in that case a pressure gradient is present at hydrostatic equilibrium. Imagine, for example, a self-gravitating gas sphere in hydrostatic equilibrium in which there is a negative radial temperature gradient, with the gas temperature decreasing radially outward from the center of the sphere. An application of the SPH equations, with the AC term of Eq. (22) now using the signal velocity (24), leads to a heat flux from the inner to the outer regions and, in the long term, to the development of an unphysical isothermal temperature profile. A signal velocity that avoids this difficulty is given by (25)which corresponds to the formulation proposed by Wadsley et al. (2008) in which a dissipative term is added to the evolution of the thermal equation with the purpose of modeling heat diffusion due to turbulence. The new term is constructed in analogy with the subgrid-scale model of Smagorinsky (1963) and is given by (26)where C is a coefficient of order unity, whose precise value depends on the problem under consideration. For the rising hot bubble problem considered by Wadsley et al. (2008), the closest agreement with the corresponding PPM results taken as reference is obtained by setting C = 0.1, higher values being too diffusive. Throughout this paper, the SPH simulations of the hydrodynamic tests are performed by adopting the expression (18) for the thermal energy dissipative term. With respect to the formulation of Wadsley et al. (2008), this approach has the advantage of using a diffusion parameter that is not constant but evolves in time according to a source term, so that the amount of diffusion away from discontinuities is minimized. For the AC signal velocity, we then use Eq. (25), whose performances in the AC formulation (18) has not been fully tested before in SPH simulations of hydrodynamic test problems.

We note that, using the signal velocity (25), the Von Neumann stability criterion becomes unimportant with respect the Courant condition (11). The Von Neumann constraint requires Δt ≤ 0.5Δx2/D, which in SPH is given by (27)Since αC ≃ O(1), this condition implies that in both the supersonic and subsonic regimes.

3. Hydrodynamic tests

In the following, simulation results obtained by applying the new AC-SPH code to a number of hydrodynamic test problems are discussed with the objective of validating the code and assessing its performance. To this end, the simulation results of the tests will be compared with the corresponding ones obtained by previous authors using different codes and/or numerical methods. The problems considered are usually presented in the literature in dimensional or complexity order, but here we follow a different approach. Since the KH instability is the hydrodynamic test that initially (Agertz et al. 2007) inspired the discussion about the difficulties of standard SPH in properly handling the development of instabilities, and it is also the most widely considered in this context (Abell 2011; Price 2008; Cha et al. 2010; Heß & Springel 2010; Read et al. 2010; Valcke et al. 2010; Murante et al. 2011), we here discuss first the two-dimensional (2D) version of the test in detail. The setup of the other runs is then considered in the light of the results obtained from the 2D KH test.

3.1. 2D Kelvin-Helmholtz instability

The KH instability arises in the presence of shear flow between two fluid layers when a small velocity perturbation is imposed in the direction perpendicular to the interface between the two fluids. The development of the instability is characterized by an initial phase, where the fluids interpenetrate each other, and then the forming of vortices, which become progressively more pronounced and turn into KH rolls at the onset of non-linearity. In the case of incompressible fluids for a sinusoidal perturbation of wavelength λ, this phase is reached on a timescale (Chandrasekhar 1961) (28)where ρ1 and ρ2 are the fluid densities and v = v1 − v2 is the relative shear velocity. A proper modeling of instabilities in numerical hydrodynamics is essential since the KH instability plays a crucial role in the development of turbulence that occurs in many hydrodynamical phenomena. In particular, the KH instability is relevant to many astrophysical processes, such as gas stripping from cold gas clouds occurring in galactic halos and the production of entropy in the intracluster gas of galaxy clusters due to the injection of turbulence.

The growth of the KH instability in hydrodynamic simulations has been addressed by various authors (Abell 2011; Price 2008; Wadsley et al. 2008; Read et al. 2010; Valcke et al. 2010; Heß & Springel 2010; Cha et al. 2010; Murante et al. 2011). These studies found that the development of the instability is artificially suppressed because of two distinct effects: the first problem is the so-called local mixing instability (LMI) due to entropy conservation, which inhibits mixing on the kernel scale and thus introduces pressure discontinuities; the second problem is caused by errors in the momentum equation, which cannot be reduced by increasing the number of neighbor particles without causing particle clumping.

Given the wide variety of numerical parameters and initial conditions with which the KH instability has been addressed, we choose here to perform the tests using as reference the simulations presented by Valcke et al. (2010). In particular, the authors implement in their SPH equations an AC term similar to that of Eq. (18), but with the AC parameter set to unity, and with a signal velocity given by Eq. (24). The authors performed a systematic analysis of the capabilities of SPH to capture the KH instability using different SPH formulations, kernels, numerical resolutions, and KH timescales. This choice allows us to assess the code capabilities by constructing a large suite of simulations whose numerical results can be contrasted against the corresponding ones discussed by Valcke et al. (2010).

3.1.1. Initial conditions set-up

The problem domain consists of a periodic box of unit length with Cartesian coordinates x ∈  { 0,1 } , y ∈  { 0,1 } . Within this domain, there is a fluid with adiabatic index γ = 5/3, which satisfies the conditions (29)As in Valcke et al. (2010), we choose here ρ1 = 10 and ρ2 = 1, with the index 1 referring to the high density layer. This choice is motivated by the finding that the difficulties of SPH in reproducing KH instabilities increase as the density contrast between the two fluid layers gets higher. The two layers are in pressure equilibrium with P1 = P2 = 10, so that the sound velocities in the two layers are and , respectively. The layers slide against each other with opposite shearing velocities v1 =  −v2 and to seed the KH instability a small single-mode velocity perturbation is imposed along the y-direction (30)where w0 = 0.025 and λ = 1/6. To restrict the perturbation to spatial regions in the proximity of the interfaces, the perturbation (30) is applied only if  | y − σ |  < 0.025, where σ takes the values of 0.25 and 0.75, respectively. With this choice of parameters, the Mach number of the high-density layer is M ≃ v1/c1 ≃ v1/1.29 and the KH timescale is τKH ≃ 0.29/v1. As in Valcke et al. (2010), we consider KH simulations with a range of five different Mach numbers; Table 1 reports the values of M together with the respective values of v1 and τKH.

Table 1

KH parameters for the simulations.

To implement the density set-up (29), a two-dimensional lattice of equal-mass particles is placed inside the simulation box. We adopt here an isotropic hexagonal-close-packed (HCP) configuration for the particle coordinates instead of the more commonly employed Cartesian grid. The advantage of this configuration is that, for a given number of neighbors, it gives a more robust density estimate owing to its symmetry properties. To construct the initial density configuration (29), the lattice spacing of the particles is varied until the SPH density estimate (1) satisfies the required values within a certain tolerance criterion ( ≲ 1%) for the relative density error. The simulations were run using a total number of N = 5122 particles and the initial particle specific energies were assigned after the density calculation so as to satisfy pressure equilibrium. This particle number is larger than that used in the runs of Valcke et al. (2010) ( ≃ 200 K), but guarantees a density setup with the specified tolerance criterion.

As noticed by Valcke et al. (2010) and other authors, SPH is a numerical method that can only represent smoothed quantities, hence applying it to hydrodynamic problems where strong density gradients are present can lead to inconsistencies This is, in fact, the situation encountered by standard SPH in the treatment of KH instabilities, where the lack of entropy mixing induces an artificial pressure discontinuity at fluid interfaces with a jump in density.

Motivated by these difficulties, we consider here a set of runs in which the density discontinuity at the interfaces is replaced by a smooth transition. To allow us to make a consistent comparison with the corresponding runs of Valcke et al. (2010), we adopt the same smoothing profile (31)where the coefficients are given by and y′ = y − σ + δ/2. In Eq. (31), the sign in front of the coefficient A refers to σ = 0.25, 0.75, respectively. The parameters β and δ determine the thickness of the density transition and take the values β = 10, δ = 0.5/0.75.

The following procedure is adopted to construct a lattice configuration that satisfies the density profile given by Eq. (31). An HCP lattice of particle is first constructed in the domain 0 ≤ x < 1 and 0 ≤ y ≤ 0.5 with a spacing adown such that the SPH density is ρ = ρ2. Proceeding upward from y = 0, the lattice spacing is progressively reduced according to a = adown [ρ2/ρ(y)] 1/2 until y = 0.25 and Ndown particles are left. The same procedure is applied to a high-density lattice that in the same domain satisfies the condition ρ = ρ1: starting from y = 0.5 and proceeding downward, the lattice spacing is increased so that a = aup [ρ1/ρ(y)] 1/2 until y = 0.25 and Nup particles remain of the original high-density lattice. The whole procedure is numerically iterated until the numbers Ndown and Nup satisfy the conditions Ndown + Nup = N/2 and Nup/Ndown = ρ1/ρ2. The lattice is then replicated in the top half of the box (y → 1 − y) so that the initial conditions are fully symmetric around y = 0.5. In the following, simulations with initial conditions for which the particle positions are arranged in a unsmoothed HCP lattice are denoted with the suffix SPH, whereas all of the others adopt the smoothed density profile constructed according to the procedure described here.

Another issue that is found to have a significant impact on the ability of standard SPH to properly capture KH instabilities is the choice of the kernel function. According to Read et al. (2010), the accuracy of the momentum equation for particle i is governed by the leading error (36)which vanishes in the continuum limit. However, this does not hold for a finite number of irregularly distributed particles or, more specifically, in the proximity of a contact discontinuity where a density step is present. A possible solution consists of increasing the number of neighbors so as to improve kernel sampling, although this approach presents some difficulties when the commonly employed M4 or cubic spline (CS) SPH kernel is used. This occurs because for a large number of neighbors the CS kernel is subject to the so-called clumping instability, in which pairs of particles with interparticle distance q = r/h < 2/3 remain close together because the CS kernel gradient tends to zero below this threshold distance. A stability analysis (Morris 1996; Børve et al. 2004; Price 2005a; Read et al. 2010) shows that for the CS kernel, a Cartesian lattice of particles is unstable for η ≃ 1.5 or Nsph ~ 28,110 when D = 2 and D = 3, respectively. The clumping degrades the spatial resolution because it reduces the effective number of neighbors with which integrals are sampled, thus still has large E0 errors even when the resolution is increased. To overcome this problem, one can modify the kernel shape in order to have a non-zero kernel derivative at the origin. However, for a fixed number of neighbors, kernels of this kind have the drawback of providing a less accurate density estimate than that given by the CS kernel since near the origin the kernels have a steeper profile.

As an alternative to the CS kernel, we consider here the linear quartic (LIQ) kernel, introduced by Valcke et al. (2010), which is a quartic polynomial for qs ≤ q ≤ ζ = 2 and is linear for 0 ≤ q < qs. The choice of parameter qs determines the quality of density estimates. From a set of 2D Sod shock tube runs, Valcke et al. (2010) recommend qs = 0.6, which is the value adopted here. For the functional form and normalization of the kernel, we refer to Valcke et al. (2010).

Another kernel that has been introduced for the purpose of avoiding particle clumping is the core triangle (CRT) kernel (Read et al. 2010), which has a constant first derivative for 0 ≤ q < α and is similar to the CS kernel for α ≤ q ≤ ζ = 2. This kernel is of the form (37)where for D = 2 and D = 3, respectively. The value of α is fixed by the condition of continuity for the second derivative, giving α = 2/3. For the grid of initial conditions previously described, the sample of KH simulations is then constructed by performing SPH runs with the same initial conditions but using different kernels. We consider the CS kernel, together with the LIQ and CRT kernels. For all of these runs, the number of neighbors is Nsph = 32   (η ≃ 1.5).

Table 2

Main simulation parameters of the KH tests.

thumbnail Fig. 1

Density maps for some of the 2D KH instability tests described in Sect. 3.1. From left to right, each column shows the panels for runs having initial conditions with the same Mach number at the corresponding time t = τKH, as listed in Table 1. From top to bottom, different rows refer to the simulations SPH, RHO, LIQ, and LP. The last one uses the linear quartic kernel but with the signal velocity (24), whereas the first three use the expression (25) (see Table 2). The plots can be compared directly with Fig. 10 of Valcke et al. (2010).

Open with DEXTER

thumbnail Fig. 2

As in Fig. 1 but for the set of simulations LIQ, CRT, and M5.

Open with DEXTER

thumbnail Fig. 3

Rendered plots of the αC parameters are shown for the same runs as in Fig. 2. The distribution of particle values  has been interpolated at the map grid points according to the SPH prescription.

Open with DEXTER

thumbnail Fig. 4

As in Figs. 1 and 2 but here at the time t = 2 for all of the panels. Note from Table 1 that this implies t ≫ τKH for runs with high Mach numbers.

Open with DEXTER

However, another solution for reducing sampling errors consists of keeping a B-spline kernel function but increasing its order (Price 2012a, Sect. 5.4). This approach has the advantage of retaining the bell shape of the kernel, a feature that provides good density estimates (Fulk & Quinn 1996). After the M4 (CS) kernel, at the next order is the M5 or quartic kernel (38)which is truncated to zero for ζ ≥ 2.5 and has σ = 96/1199π,1/20π for D = 2 and D = 3, respectively.

An additional set of SPH simulations was then performed using the M5 spline as the kernel. For consistency with the other runs, we kept the same ratio of smoothing length to particle spacing (η ~ 1.5), so that for these runs the chosen number of neighbors is Nsph = 50. A non-trivial issue concerns the role of pairing instability for this class of kernels. Because the gradient of the M5 kernel still goes to zero as q → 0, one would expect the instability still to occur for η ~ 1.5. Nonetheless, it will be seen that KH simulations with the M5 kernel do not exhibit particle clumping, in contrast to corresponding simulations performed with the CS kernel. This suggests that the stability properties of the M5 kernel are better than those of the CS one; this topic is discussed in a subsequent dedicated section (Sect. 3.1.3).

Table 2 summarizes the main simulation parameters used in the KH tests. For comparison purposes, a set of mirror runs was carried out for the LIQ simulations in which Eq. (25), for the signal velocity, was replaced in Eq. (18) by Eq. (24), which is based on pressure discontinuities.

3.1.2. Results

Figure 1 shows the density maps at t = τKH for some of the KH simulations listed in Table 2. The panels can be compared directly with the corresponding ones in Fig. 10 of Valcke et al. (2010). A visual inspection reveals that the expected features of the KH runs are correctly reproduced and of a general consistency between the results produced by the two codes. In particular, the LP simulations, which use the pressure-based AC signal velocities of Eq. (24), appear to produce results almost identical to the corresponding LIQ ones. Thus, showing how, for the KH tests examined here, the use of the two signal velocities can be considered equivalent, within the spanned range of Mach numbers.

We note, however, the absence of KH rolls for the LP run with M = 0.4, which contrasts with the same simulation in Valcke et al. (2010) where the rolls have already been developed. Given the general agreement between the two codes, it is hard to ascertain the origin of the discrepancy for this specific run. However, the panels of Fig. 1 show that the AC-SPH formulation, as discussed in Sect. 4.2, is still unable to properly capture the development of KH instabilities for very subsonic shear flows. This suggests how small differences in the time integration procedure of the two codes can manifest themselves in the long-term evolution of cold flows.

In Fig. 2, the results for the same set of KH simulations of Table 1 are shown, but with the use of kernels LIQ, CRT, and M5. An important feature is now the appearance in the CRT and M5 runs of KH rolls in the M = 0.4 case. This improvement in the description of KH instabilities suggests that the integration errors that are present with the LIQ kernel are now absent or smaller in the CRT and M5 runs. However, in the M = 0.2 case, the kernels are still unable to capture the development of the KH instability. For this test case, a high-resolution run (N = 10242) using the M5 kernel (not shown here) still reveals the absence of any KH roll at t = τKH, thus showing that in SPH the problem of an accurate description of KH instabilities in the very subsonic regimes is neither a resolution issue nor due to the AC implementation. We refer to McNally et al. (2012) for a recent discussion of this topic.

To test the effectiveness of the switch (13) in ensuring a fast response of the  parameters to the presence of thermal energy discontinuities, Fig. 3 renders the plots of the AC parameters that are shown for the test runs of Fig. 2. The maps are color-coded according to the range of values of . The highest values lie in the range  ~0.7−0.8 and are confined to a narrow strip around the interface layers, with the floor value  ~0 as the background value away from the discontinuities. We note that in general, and in particular for the M = 1 test case, the maximum values of the αC for the M5 runs are below those of the other runs. This result is due to the faster growth rate of the instability, which is ensured by an improved kernel accuracy in interpolating the variables.

The long-term evolution of the KH tests is shown at t = 2 in Fig. 4. The overall features of the density plots for different runs are similar to their counterparts displayed in Fig. 11 of Valcke et al. (2010). We note the tendency for the M5 runs to show small-scale features at the layer contacts.

The aim of this section is to test the consistency of the present AC implementation by comparing results extracted from a suite of AC-SPH simulations of the 2D KH instability problem, with those of similar runs (Valcke et al. 2010). The results of the KH tests indicate a code behavior that is in accord with that expected and the simulations of Valcke et al. (2010).

thumbnail Fig. 5

Time evolution of the velocity field amplitude in the y direction, as measured by the  = 12π mode of the Fourier transform of vy, for some of the KH instability tests described in Sect. 3.1. Each panel refers to KH simulations performed with the same Mach number, the initial conditions set-up being given in Table 1. Within each panel, the different curves are for AC-SPH runs with different simulation parameters, as specified in Table 2. The black solid line indicates the expected linear-theory growth rate.

Open with DEXTER

However, Valcke et al. (2010) argue that the absence of the AC term is not the main reason for the SPH failure to develop KH instabilities; however the presence of AC is necessary for the long-term evolution of the instabilities, but this difficulty of SPH has two distinct reasons. The first is a general problem of consistency of the initial condition set-up, as SPH by definition can only deal with smoothed quantities and therefore its application to problems where sharp discontinuities are present leads to inconsistencies. This is part of the more general problem (Robertson et al. 2010) of properly smoothing in numerical simulations of hydrodynamic test problems the discontinuities initially present at the interfaces, in order to achieve convergence in the solution. This aspect of the KH test problem can be cured by properly stretching the initial particle lattice so as to introduce a smooth density transition at the interfaces. Results suggest that there is a significant improvement in the capability of SPH to properly capture the correct growth rate of the KH instability.

The other problem that causes the poor performance of SPH when handling KH instabilities is the leading error in the momentum equation due to the incomplete kernel sampling. This error can be reduced by removing particle clumping, which depends on the kernel stability properties. In fact, the kernels LIQ and CRT were introduced (Read et al. 2010; Valcke et al. 2010) with the aim of removing the clumping instability. Since the results of the simulations suggest performances for the M5 kernel that are quite similar to those achieved by these kernels, it is therefore interesting quantify the relative performances of these kernels in a better way.

Following Junk et al. (2010), we then measure the growth rate of the KH instability for some of the runs and compare it with the linear theory growth rate expectation  ∝ et/τKH. This is achieved, using Fourier transforms, by measuring the time-evolution of the λ = 1/6 growing mode of the vy velocity perturbation component. The details of the procedure are given in Junk et al. (2010) and so are not reported here. For the sake of clarity, the results of the LIQ runs are not displayed since they are quite similar to those of the CRT simulations. Moreover, we only display the growth rates for three different KH test cases, those with Mach number M = 0.2, 0.4, and M = 1, the results of the others being intermediate between these.

A number of distinct features are apparent in the panels of Fig. 5. The first is that simulations with smoothed IC (RHO) perform systematically better than the unsmoothed ones (SPH). The second is that only for M = 1 are the simulations with kernels CRT and M5 able to correctly recover the expected growth rate. Finally, this capability progressively declines as one considers lower Mach numbers. While this behavior is in accordance with the visual impression of the maps previously displayed, and confirms that the present SPH implementation still has problems in describing instabilities in subsonic flows, we note that the performances of the M5 runs are similar to those obtained using the CRT kernel. This behavior is particularly interesting, since the simulations were performed by setting the ratio of the smoothing lengths to particle spacing to η ≃ 1.5 so that the pairing instability, which is present in runs that use the CS kernel, should be present in the M5 simulations as well.

How the clumping instability affects sampling errors can be assessed by computing the particle errors , according to Eq. (36). The distribution, throughout the simulation domain of the  errors versus y, is shown in the top panels of Fig. 6 for M = 1 simulations performed using different kernels. The solid lines represent the mean of the binned distributions. Similar plots were produced by Valcke et al. (2010) and their Fig. 3 can be used for comparative purposes.

As expected, the largest E0 errors are present in the SPH simulation, whereas better results are obtained by using the LIQ and CRT kernels, as indicated by the error distribution in their respective panels. This is a result of the absence of particle clumping for these simulations, owing to the specific stability properties of these kernels (Read et al. 2010; Valcke et al. 2010).

A striking feature is given by the E0 error distribution of the simulation performed using the M5 kernel, which shows how the magnitude of the E0 errors are even smaller than those of the two runs CRT and LIQ for this kernel. This result is in accordance with what has been found previously by analyzing the growth rates, suggesting that, since all of the simulations were performed using the same η, the stability properties of the M5 kernel are better than those of the CS one.

thumbnail Fig. 6

Top panels: distribution at time t = τKH of the errors  | Ei0 |  plotted versus y, as defined by Eq. (36), for the KH runs with Mach number M = 1 and different simulation parameters (See Table 2). The red histograms show the mean values of the binned distributions. Bottom panels: each panel shows the nearest neighbor map of the run in the corresponding top panel. This is defined as the distribution, interpolated at the map grid points of the normalized distances , where  is the distance  | xi − xk |  of the kth neighbor of the particle i and the neighbors are sorted so that .

Open with DEXTER

To further investigate this point, the bottom panels of Fig. 6 show the “nearest neighbor” map of the corresponding top panels. This is defined by interpolating the quantity  at the grid points, according to the SPH prescription, where for the k − th neighbor and the neighbors have been sorted according to their distance from the particle i itself.

The map of the M5 kernel indicates a distribution of the second nearest-neighbor distribution  that is quite similar to that displayed by the CRT and LIQ kernels. The only difference occurs at the layer interfaces, where the distribution of the quantities  for the M5 kernel is slightly shifted towards smaller values than for the other kernels. Conversely, we note that because of the pairing instability, for the CS kernel the distribution of the neighboring distances is flipped with respect to that of the other kernels.

The results of Fig. 6, however, clearly show the absence of clumping instability for the M5 kernel. In the next section, we investigate, in more detail, the stability properties of this kernel.

3.1.3. Stability issues

The stability properties of SPH have been investigated by a number of authors (Swegle & Hicks 1995; Morris 1996; Monaghan 2000; Børve et al. 2004; Read et al. 2010; Dehnen & Aly 2012). The instabilities are studied analytically by analyzing the dispersion relation for sound waves of small amplitude, propagating in a uniform medium. We now derive the dispersion relation for the SPH equation of motion in a manner similar to the analysis of previous authors. A uniform lattice of equal-mass particles of mass m, density ρ0, pressure P0, and sound speed is perturbed by a small wave (39)where a is the perturbation, k is the wavevector, and  are the unperturbed particle positions. By linearizing the equation of motion for the perturbation, one has the dispersion relation where the summations are over the neighbors j of particle i, the vector qi is defined as (40)and H(W) is the Hessian of the kernel (41)For a given smoothing length h and wavevector k, the stable solutions of Eq. (40) are defined by the condition ω2 ≥ 0. Solutions for which ω2 < 0 represent exponentially growing or decaying perturbations. Moreover, it is useful to define a numerical sound speed and a scaled numerical sound speed . The condition should clearly be satisfied for the sound speed to be correctly modeled.

thumbnail Fig. 7

Contour plots of the dispersion relation given in Eq. (40) are shown as a function of wavenumber k and smoothing length h for a regular lattice of particles with spacing Δ, where we consider a wavevector with orientation k = k(1,0,0). Gray areas indicate the instability regions for which ω2 < 0. From the left to right, the top panels show in the case of the M5 kernel the instability regions of the longitudinal and the two transverse waves. The bottom panels refer to the CS (M4) kernel.

Open with DEXTER

We now examine the stability properties of the CS and M5 kernels. For simplicity, we consider a plane wave propagating along the x-axis, k = k(1,0,0), and assume γ = 5/3. The bottom panels of Fig. 7 show, for the longitudinal and the two transverse waves of the perturbation, the instability regions of the CS kernel, which are denoted by the gray areas and represent the solutions to Eq. (40) in the domain (h,k), for which ω2 < 0. Similarly, the top panels are for the M5 kernel.

The longitudinal wave perturbation is responsible for the clumping instability, whereas the traverse waves give rise to the so-called banding instability (Read et al. 2010). Unlike the clumping instability, the latter is relatively unimportant in causing sampling errors (Read et al. 2010) and is not considered further here. A comparison of the two stability plots for the longitudinal wave solution clearly shows that the stability properties of the M5 kernel are much better than those of the CS one.

This is part of a more general result that was previously recognized by Morris (1996): the stability properties of SPH improve when higher-order spline kernels are used in place of the CS kernel. This occurs, basically, because the higher the order of the spline, the better it approximates a Gaussian kernel. In Eq. (40), one can see that the numerical sound speed  depends on the first and second derivatives of kernel W. Ideally, one should have to accurately describe the sound wave propagation and this condition is fulfilled when smoother kernels are used, so as to keep the numerical dependence of  as weak as possible.

We note however that here, in contrast to the stability properties of the CRT kernel, the clumping instability is not entirely suppressed but instead it is present whenever η ≳ 2.5.

Finally, it must be stressed that the paring instability that occurs in the hydrodynamic tests described here, is unlikely to have a significant impact on many astrophysical problems of interest, where very cold flows are generally absent. Rasio & Lombardi (1999) estimate for instance, from SPH simulations of a stationary fluid, that lattice effects become important when velocity dispersions are below  ~3%−4% of the sound speed.

The results of this section and the previous one, therefore, suggest that to avoid the pairing instability, the M5 kernel can be considered a viable alternative to the use of the otherwise steeper CRT and LIQ kernels, provided that the parameter η is consistently rescaled. In the next section, we then discuss the related performances of these kernels in a test simulation.

3.2. Sod’s shock tube

A classic test used to investigate the hydrodynamic capabilities of SPH codes is the Sod shock-tube problem (Hernquist & Katz 1989; Wadsley et al. 2004; Springel 2005; Price 2008; Tasker et al. 2008; Rosswog 2009). This test consists in a fluid, initially at rest, in which a membrane located at x = 0 separates the fluid on the right, of high density and pressure, from the fluid on the left, of relatively lower density and pressure. The membrane is removed at t = 0 and a shock wave develops propagating toward the left, followed by a contact discontinuity and a rarefaction wave propagating to the right.

A well-known problem with standard SPH codes in reproducing the analytic solution of the shock-tube problem is the presence of a pressure discontinuity that arises across the propagating contact discontinuity. Simulations incorporating an artificial thermal conductivity term in the SPH equations (Price 2008; Rosswog 2009; Price 2012a) show shock profiles in which density and thermal energy are resolved across the discontinuity, hence producing a continuous pressure profile. However, in these runs the AC formulation adopts the AC signal velocity (24) where jumps in thermal energy are smoothed in accordance with the presence of pressure discontinuities. This is in contrast to the AC signal velocity (25) employed here, for which in the absence of shear flows, contact discontinuities are unaffected and therefore cannot be used to remove the blip seen at the contact discontinuity in the pressure profile of the shock tube SPH runs.

The application of the AC-SPH code to this test is nonetheless of interest because it can still be used to validate the code’s performances. In particular, we consider a 3D setup of the shock-tube test and with these initial conditions we construct a set of AC-SPH simulations performed with different kernels and neighbor numbers. Shock tube profiles extracted from these simulations are compared with the aim of assessing the goodness of different kernels in reproducing, for a given test problem, profiles of known analytic solutions.

The initial condition setup consists of an ideal fluid with γ = 5/3, initially at rest at t = 0. An interface set at the origin separates the fluid on the right of density and pressure (ρ1,P1) = (4,1) from the fluid on the left with (ρ2,P2) = (1,0.1795). To construct these initial conditions, a cubic box of side-length unity was filled with 106 equal mass particles, so that 800   000 were placed in the right-half of the cube and 200   000 in the left-half. The particles were extracted from two independent, uniform glass-like distributions contained in a unit box consisting of 1.6    ×    106 and 4 × 105 particles, respectively. This initial condition setup is the same as that previously implemented in Sect. 5.1 of V11, to test the time-dependent AV scheme described in Sect. 2.2.

thumbnail Fig. 8

Results at t = 0.12 of the 3D shock tube test. The profiles of density, pressure, thermal energy, and velocity are plotted clockwise from top left. The solid black line represents the analytical solution, while lines with different styles and colors are the profiles of the AC-SPH runs with different kernels and neighbor numbers, as illustrated in the pressure panel.

Open with DEXTER

For the same initial setup, to investigate the performances of different kernels in reproducing the analytic profiles of the shock-tube problem, we perform runs with different kernels and neighbor numbers. The kernels considered are CS, LIQ, CRT, and M5. The number of neighbors varies as a power of two between 32 and 128. The SPH runs were realized by imposing periodic boundary conditions along the y and z axes and the results were examined at t = 0.12. We show the results obtained using the standard AV formulation, the results of the other runs being unimportant from the viewpoint of the kernel performances.

Some of the simulation profiles extracted from the 3D SPH runs of the shock-tube test are shown in Fig. 8. A striking feature is the large scatter between the pressure profiles of the runs performed using different kernels or neighbor numbers. The same behavior is present for the thermal energy profiles, whilst very similar density profiles are exhibited by the same runs.

There are several conclusions that can be drawn from the results of Fig. 8. The performances of the M5-64 run are quite similar to those of CS-32, although for the former simulation a closer inspection reveals a slightly better treatment of the thermal energy spike at the contact discontinuity, the spike being due to the initial condition set-up.

The worst results are obtained by the LIQ simulations when using Nsph = 32 or Nsph = 64 neighbors. The profiles of the LIQ-128 run are not shown here to avoid overcrowding in the plots, and are quite similar to those of the simulation CRT-128. This clearly demonstrates the need for this class of kernels to use a large (say  ≳ 128) number of neighbors, so as to compensate for the density underestimate due to the steeper profiles introduced to avoid the pairing instability.

To better quantify this density bias, Table 3 reports, for different kernels and neighbor numbers, the mean SPH density estimated from a glass-like configuration of N = 106 particles, in a unit periodic box of total mass one. The results illustrate how the density estimate of the M5 kernel with 64 neighbors is comparable to the CS one obtained using 32 neighbors, the value of η being the same (η ≃ 1). Conversely, for the LIQ and CRT kernels, only when Nsph = 128 does the mean density approach the M5-64 estimate.

Table 3

Average densities and sample standard deviations estimated from the SPH densities of a configuration of one million particles.

The density values given in Table 3 also help us to explain the rather poor performances of the LIQ kernel when using 32 neighbors. From Fig. 8, it can be seen that for the corresponding run the relative pressure error is εP ≃ 10% in the unperturbed right zone of the cube. This error is already present in the pressure profile at t = 0 and originates from the density error due to the adopted kernel and neighbor number, together with the use of an entropy-conserving code to perform the simulations. Particle entropies are initially assigned according to the initial conditions so that, for particles satisfying xi > 0, . During the integration, particle pressures are calculated according to , and for x > 0 a relative pressure error εP ≃ γερ ≃ 10% is present when ερ ≃ 6%.

Taken at face value, the results of Table 3 demonstrate that in 3D simulations a conservative lower limit for the kernels with a modified shape should be to assume at least Nsph ≳ 200 neighbors. In fact, Read & Hayfield (2012) present a new formulation of SPH where they adopted, as reference, the so-called high-order core triangle (HOCT, Price 2005a). The profile of this kernel is a generalization of that of the CRT kernel and the authors assumed Nsph = 442 as the reference value for their scheme.

However, the results of the previous sections suggest that the stability properties of the M5 kernel can be profitably used, with an appropriate choice of the η parameter, to avoid the pairing instability when dealing with tests of hydrodynamic instabilities in which cold flows are present.

Finally, the density estimates of Table 3 suggest that great care should be taken when deciding the goodness of a particular kernel based on its relative performance in terms of density estimates. The results of the 3D Sod shock tube clearly indicate how there could be a large difference between the simulation and the expected solution profile of some hydrodynamic variable, such as pressure or thermal energy, and, at the same time, a much smaller difference in the corresponding density profile.

3.3. Sedov blast wave

The Sedov blast-wave test is used to validate, in three dimensions, the code capability in the strong shock regime. The test consists of a certain amount of energy E being injected at t = 0 into a very small volume of an ambient medium of uniform density ρ. The spherically symmetric shock propagates outward from the initial volume and at time t the shock front is located at the radius (Sedov 1959) (42)where β ~ 1.15 for γ = 5/3.

Previous investigations (Rosswog & Price 2007; Merlin et al. 2010) showed that, owing to the large discontinuities initially present in the thermal energy, incorporating an artificial conduction term in the energy equation greatly improves the description of the shock front in the simulations. Without the presence of this term, the initially large discontinuity in the thermal energy soon gives rise to a disordered particle distribution that degrades the shock profile (Rosswog & Price 2007).

The initial setting of the test was realized as follows. A HCP lattice of N = 2 × 643 equal mass particles was arranged in a cubic box of side length unity. The particle masses were chosen so as to give ρ = 1 and periodic boundary conditions were imposed. The nearest particle to the position  { 0.5,0.5,0.5 }  was chosen as the center particle. To consistently represent a point-like explosion with the given numerical resolution, the particles j comprised within the kernel radius ζhi of the center particle i were given an initial thermal energy such that the total injected energy was E = 1. This blast wave energy was distributed among the neighboring particles not uniformly but with a weight proportional to Wij.

Using this initial condition set-up, we ran three test simulations, which differed in the choice of the adopted kernel and neighbor number. For the CS kernel, we used Nsph = 64 neighbors. We also ran two other test cases, now using the CRT and M5 kernel and Nsph = 128 neighbors. All of the simulations were performed using the implemented AC scheme and the time-dependent AV formulation with parameters given by  {αmin,αmax,ld}  =  {0.1,1.5,0.2}.

Figure 9 shows the radial density profiles at t = 0.06 for the different test runs. The solid black line represents the expected analytic solution, with the shock front being located at r ≃ 0.37. Radial simulation profiles were obtained by averaging, for each radial bin rk, SPH densities calculated from the particle distributions over a set of (θ,φ) = (20,20) grid points uniformly spaced in angular coordinates: these were located at the surface of a sphere with radius rk. The radial spacing was not uniform but was chosen so as to guarantee an accurate sampling of density in the proximity of the shock front, with about  ~40 radial bins between r ~ 0.34 and r ~ 0.42.

All of the simulations are in fair agreement with the analytical solution and the differences between the simulation profiles are negligible. At the shock front, the profiles exhibit a density jump of about  ~2, whereas the analytical solution gives a compression factor of γ + 1/γ − 1 = 4. These results are in accordance with previous findings (Rosswog & Price 2007; Springel 2010b; Heß & Springel 2010) and indicate that for 3D SPH simulations of the Sedov-Taylor point explosion problem the simulation profiles converge to the analytical solution as the resolution is increased, with approximately  ~3453 particles (Rosswog & Price 2007) being required to fully resolve the shock front.

Figure 9 can also be used to verify the behavior of the individual timestep algorithm. For problems involving very strong shocks, as demonstrated by Saitoh & Makino (2009), individual particle timesteps must be properly restricted so as to avoid in the proximity of the shock front their failure to satisfy the local Courant condition, thus leading to inaccuracies in the integration. The SPH simulations using individual timesteps, but without an appropriate limiter, fail to predict the expected solution profile (see Fig. 3 of Saitoh & Makino 2009). The profiles of Fig. 9 demonstrate the correctness of the algorithm used to update the timesteps, although different from the one devised by Saitoh & Makino (2009).

The present scheme adopts individual particle timesteps Δti that can vary in power-of-two subdivisions of the largest allowed timestep Δt0 ≥ Δti (Hernquist & Katz 1989). At each step, particles whose time bin is synchronized with the current time are defined as “active” and their hydrodynamic quantities, as well as their smoothing lengths and densities are consistently updated. However, non-active particles j that are neighbors of an active particle i are defined as “low-order active” particles, and their hydrodynamic variables, as well as their time-step constraints but not their forces, are recalculated. This criterion is applied regardless of the local shock conditions and no particular conditions are imposed on Δtj that depend on Δti.

thumbnail Fig. 9

Radial density profiles of the 3D Sedov blast wave test at t = 0.06. The solid black line indicates the analytic solution, while the simulation profiles are obtained by averaging for each radial bin SPH densities calculated from the particle distributions over a set of grid points located at the surface of a spherical shell and uniformly spaced in angular coordinates. All of the runs were performed using N = 524   288 particles.

Open with DEXTER

thumbnail Fig. 10

Density maps of the blob test in the central plane x = Lx/2 at t = 1, 2, 3 for SPH runs with (CS, M5) and without (CS NOAC) the AC term. Time is in units of τKH ~ 1.6τcr and the number accompanying the kernel label indicates the number of neighbors of the run. Axis units are in Mpc.

Open with DEXTER

3.4. The blob test

The blob test is another hydrodynamic test where the results of standard SPH differ significantly from those produced by grid-based simulations (Agertz et al. 2007; Read et al. 2010; Cha et al. 2010; Heß & Springel 2010; Murante et al. 2011). The test consists of a gas cloud of radius Rcl placed in an external medium ten times hotter and less dense than the cloud, so as to satisfy the pressure equilibrium. A large enough wind velocity vW is given to the hot low-density medium so that a strong shock wave strikes the cloud. The interaction of the cloud with the supersonic medium produces a number of effects that are of interest in an astrophysical context, such as gas stripping and fragmentation.

The blob is initially perturbed by the development of Richtmyer-Meshkov and Rayleigh-Taylor instabilities (Agertz et al. 2007). Afterwards, large-scale (~Rcl) KH instabilities are created at the cloud surface because of the shear flows caused by the supersonic wind. This non-linear phase is supposed to develop over a timescale (Agertz et al. 2007), where χ is the density contrast, after which cloud disruption will take place.

To investigate the capability of the AC-SPH code to properly follow the hydrodynamics of the blob test, we compare results extracted from our set of SPH simulations realized with the same initial conditions but with different numerical parameters. The numerical setup of this test is the same as in Read et al. (2010), to which we refer for more details. A spherical cloud of radius Rcl = 197 kpc is placed in a periodic, rectangular box of size  { Lx,Ly,Lz }  =  { 2,2,8 }  Mpc. The cloud has density ρcl = 4.74 × 10-33   gr   cm-3 = χρext and temperature Tcl = 106   K = Text/χ. The ambient medium has density and temperature so that χ = 10. The cloud is initially located at  { 1,1,1 }  Mpc and the ambient medium is given a wind velocity vW = 1000   km   s-1, so that for an adiabatic index γ = 5/3 its Mach number is M = 2.7. An HCP lattice of equal mass particles is constructed to satisfy the above density requirements so as to use for this version of the test a total number of N ~ 1.1 × 106 particles, as in Heß & Springel (2010). Finally, a velocity perturbation is imposed at the cloud surface in order to trigger the development of an instability; the amplitude and modes are given in Appendix B of Read et al. (2010).

We compare the results of SPH simulations where three different spline kernels were used CS, CRT, and M5. We use 128 neighbors for the simulations with the CS and CRT kernels and 220 for the run with the M5 kernel, so that the ratio η is the same for all the runs (η ~ 1.5). For the CS kernel, we run a simulation where the AC term given by Eq. (18) is absent in the equation of thermal energy evolution (CS-128 NOAC), and three other simulations (CS-128, CRT-128, M5-220) in which the AC term (18) is incorporated in the energy evolution equation. In the following, the simulation AV parameters are set as follows:  { αmin,αmax,ld }  =  { 0.1,1.5,0.2 } .

For several runs, Fig. 10 shows the gas density maps at three different times t = 1, 2, 3, the time being in units of τKH ~ 1.6τcr (Agertz et al. 2007). The maps have been evaluated according to the SPH prescription on a 2D grid of (200x800) points in the central yz plane located at x = Lx/2. The top panels of Fig. 10 are for the standard NOAC run CS-128. It can be seen that in this case, unlike the results of mesh-based simulations (Agertz et al. 2007), the cloud survives the impact of the supersonic wind striking its surface and there is no disruption. This is due to the absence of fluid instabilities developing at the cloud surface. If they were present, they would in turn lead to the stripping of material and the cloud break-up. This code behavior is similar to what was found in previous findings (Agertz et al. 2007; Read et al. 2010; Heß & Springel 2010).

thumbnail Fig. 11

Mass loss of the cloud versus time for the blob test simulations. Particles are defined as a cloud member if their temperatures and densities fulfill the conditions Ti < 0.9Text and ρi > 0.64ρcl. Left panel: simulations have been performed incorporating the AC term into the SPH equations and using the kernels CS, CRT, and M5. For the CS kernel, a simulation is also shown as reference without AC (NOAC, solid black line). Right panel: for the CS kernel, simulations with additional constraints on the particle AC parameters  have been performed using N = 106 (N6) and N = 105 (N5) particles, see text for more details. Black squares have been extracted from Fig. 10 of Springel (2011) and indicate the behavior of the cloud mass versus time in a similar test performed using the moving-mesh code AREPO and with (64 × 64 × 128) resolution elements.

Open with DEXTER

In contrast, incorporating the AC diffusion term in the SPH equations leads to a significant improvement in the code capability to properly model the cloud evolution. This can be seen from the middle and bottom row panels of Fig. 10, in which the density maps are displayed for the CS-128 and M5-220 runs. The simulation performed using the CRT kernel is not shown since its performances are quite similar to those obtained using the CS kernel.

This behavior is consistent with the expectation that introducing AC removes the numerical effects that in SPH prevent the treatment of contact discontinuities when large density jumps are present, thus the inconsistencies that suppress the growth of the instabilities. The AC-SPH formulation presented here can therefore, at least qualitatively, correctly follow the time evolution of the cloud as in the other SPH schemes that have been proposed (Read et al. 2010; Heß & Springel 2010; Murante et al. 2011; Saitoh & Makino 2012).

To quantify the code performances in a more quantitative way, Fig. 11 (left panel) shows the mass loss of the cloud as a function of time for the different runs. We follow previous definitions (Agertz et al. 2007) and a gas particle is defined to still be a member of the cloud if at the considered epoch its gas and temperature fulfill the conditions ρ > 0.64ρcl and T < 0.9Text. The standard SPH results show a cloud that is not disrupted and whose mass, after an initial transient, stays constant at about half of the original value. A striking result is instead given by the runs employing the AC term. For these simulations, there is now a high mass-loss rate occurring at early times, followed by complete cloud disruption at t ≳ 3τKH.

We note that the mass loss rate of the AC-SPH runs does not depend in a significant way on the choice of the kernel. This suggests that the cloud disruption is driven by large-scale instabilities and is relatively insensitive to small-scale perturbations. Given the similarities displayed in Fig. 11, left panel, by the mass loss rates of the AC-SPH simulations employing different kernels, when referring to the left panel we adopt the term mass loss rate to generically indicate the behavior of these curves.

In the test considered here, the new AC scheme is clearly now capable of properly removing the surface effects, present across the contact discontinuity in the standard SPH version, which artificially suppress the growth of hydrodynamic instabilities. However, a comparison of the mass loss rate with the corresponding one produced using the new moving-mesh code AREPO in simulations of equivalent resolution (Springel 2011, Fig. 10) shows that for the AC runs presented here the mass depletion of the cloud occurs much faster. This discrepancy suggests that the processes of heat diffusion, which in the adopted numerical scheme are mediated via the AC parameter αC, should be somehow constrained by a physically motivated mechanism, which was not considered in the discussion of Sect. 2.3. This mechanism should be introduced with the purpose of avoiding a heat transfer mechanism, as governed in the code by Eq. (18), which is overly diffusive.

To evaluate the relative effectiveness of the heat and momentum transport, in the theory of heat transfer the Prandtl number Pr is defined as the ratio of the kinematic viscosity ν to the thermal diffusion coefficient D, Pr = ν/D (Blundell & Blundell 2006). For gases, the transport coefficients for the transport of heat and momentum are nearly equal and the Prandtl number is of order unity, with when γ = 5/3 (Blundell & Blundell 2006). This suggests that a constraint on the particle AC parameter  can be obtained by setting (43)where the equivalent kinematic viscosity coefficient νAV due to AV is given approximately by (Lodato & Price 2010), and the numerical heat diffusion coefficient  has been estimated according to Eq. (21). We note that the reverse of Eq. (43) is relatively unimportant, since the source term SCi is driven by a second order derivative, while for AV it is a first order derivative that determines Si.

It can be seen from Eq. (43) that, for the AC signal velocity adopted here, the constraint on the numerical heat diffusion becomes important in the presence of supersonic flows. Moreover, the condition Pr ~ O(1) is valid only for gases, where the transport of momentum and energy occurs simultaneously. For liquids, according to the dominant mechanism of heat conduction, the value of the Prandtl number can vary by several orders of magnitude among different substances (Dimotakis 2005). The condition given by Eq. (43) should, therefore, be considered problem-dependent.

A rigorous procedure for deriving the upper limits on the set of parameters at any given timestep would require us to first define for the particle i the numerical kinematic viscosity and heat diffusion coefficients (44)where the bulk viscosity ζ is 5/3 of the shear viscosity η ≡ νρ (Lodato & Price 2010) and it is understood that the SPH expressions should be used for the operators at the denominators.

The upper limits on the parameters are then given by the conditions (45)which must be simultaneously satisfied by all the parameters. In place of the procedure just described, we adopt here a simpler approach and using Eq. (43) we estimate the upper limits on  as (46)where the notation j <  means that the maximum is taken for all the neighbors j that satisfy the condition vij·rij < 0. This is to be consistent with the definition of AV, for which Πij is non-zero only for approaching particles. Clearly, this definition does not guarantee an accurate constraint of the  parameter, but the use of the maximum among all the approaching pairs should provide a floor value for the effective constraint.

For the CS kernel, we then performed a simulation identical to the one previously discussed, but now imposing the additional constraint given by Eq. (46) on the αC parameters. This simulation is labeled as MAX N6 and the corresponding mass-loss curve is displayed in the right-hand panel of Fig. 11. To assess resolution effects, we ran a mirror simulation now using N = 105 particles (MIN N5). Moreover, the uncertainties associated with the use of Eq. (46) can be estimated by looking at the mass loss rate of a simulation (MIN N6) where instead of Eq. (46) the same equation was used but the constraint was derived using the MIN operator instead of the MAX one.

For comparative purposes, in the right-hand panel of Fig. 11 the mass loss of the cloud (black squares), as found in a similar blob test performed by Springel (2011) using the new moving-mesh code AREPO, has been also inserted. The number of resolution elements is 64 × 64 × 128, so that the resolution is approximately equivalent to that of the tests shown here.

A striking result seen from the behavior of the mass loss rates is that, introducing the constraint (46) on the numerical heat transfer in the simulations, the mass depletion of the cloud is now in better agreement with previous results and in particular with the cloud mass evolution produced in the blob test by the Voronoi-based code AREPO.

Uncertainties associated with the procedure adopted to estimate the constraint on the αC parameters are of limited impact, as can be assessed by the differences between the mass loss rates exhibited by the runs MAX N6 and MIN N6. Moreover, for the run MAX N5 the mass loss rate is much stronger than in the case MAX N6. This is indicative of the numerical diffusion effects that dominate the blob evolution. A numerical simulation with N = 107 particles is clearly required to clarify this issue; however, the agreement with previous results of the mass loss for the run MAX N6 suggests that convergence is achieved using as few as N = 106 particles.

Finally, both the runs MAX N6 and MIN N6 exhibit at t ≳ 3τKH a remnant cloud mass that is  ~20−30% of the initial mass. Although this result can still be due to the use of a constraint on the numerical heat diffusion that still needs to be refined, we note that these remaining masses are equally present in SPH blob tests that have adopted, with the purpose of avoiding the problems of standard SPH, completely different approaches (Heß & Springel 2010; Murante et al. 2011).

We therefore suggest that this underestimate in the stripping rate of the blob mass is not due to the AC scheme, but depends rather on the SPH formulation and is associated with the intrinsic errors in gradient estimates. These errors, in turn, lead to the suppression of the small-scale instabilities. Such a issue will be discussed further in detail in Sect. 4.2, which is dedicated to the Rayleigh-Taylor instability.

To summarize, SPH simulations that incorporate the AC scheme can successfully be used to accurately model the blob evolution. However, the results of the tests indicate that it is necessary to properly constrain the AC parameter αC in order to avoid unphysical heat diffusion when strong supersonic flows are present.

4. Gravity tests

To validate the performances of the AC signal velocity (25), we investigate several hydrodynamical test problems in which the self-gravity must be necessarily taken into account to properly model the system evolution. We first considered the 3D collapse of a cold gas sphere initially at rest. This test is customarily used in SPH to judge the code capability when strong shocks are present, such as those occurring during the formation of self-gravitating structures. We then examined a 2D version of the classic Rayleigh-Taylor instability test and finally the new code was used to assess the behavior of cluster entropies in cosmological structure formation.

4.1. Cold gas sphere

A standard hydrodynamical test for SPH codes in which gasdynamics is modeled including self-gravity is the 3D collapse of a cold gas sphere, also commonly referred to as the “Evrard” collapse test (Evrard 1988). The test follows in time the adiabatic collapse of a initially cold gas sphere and has been widely used by many authors (Hernquist & Katz 1989; Steinmetz & Mueller 1993; Wadsley et al. 2004; Springel 2005; Wetzstein et al. 2009; Springel 2010b; Heß & Springel 2010, V11) as a standard test for SPH codes.

thumbnail Fig. 12

Radially averaged profiles at t = 0.8 of the Evrard collapse test. Clockwise from the top left panel: profiles of density, entropy, artificial conductivity parameter αC, and radial velocity. Curves with different line styles and colors refer to SPH runs performed using different kernels and with (AC) or without (NOAC) the AC term in the SPH energy equation. The black solid lines indicate the 1D PPM reference solution of Steinmetz & Mueller (1993).

Open with DEXTER

The gas cloud is spherically symmetric and initially at rest with mass M, radius R, and density profile (47)The gas obeys the ideal gas equation of state with γ = 5/3 and the thermal energy per unit mass is initially set to u = 0.05GM/R. The SPH simulations are performed using units for which G = M = R = 1 and the chosen time unit is the cloud free-fall timescale .

With these initial conditions, the pressure support of the gas sphere is negligible and the cloud begins to collapse until a bounce occurs in the core with a subsequent shock wave propagating outward. Most of the kinetic energy is converted into heat at the epoch of maximum compression of the gas, which occurs at t ~ 1.1. The initial conditions setup is realized by stretching the radial coordinates of a glass-like uniform distribution of N = 88   000 equal mass particles located within a sphere of unit radius, so as to generate the density profile of Eq. (47).

To construct the set of SPH simulations, four tests were run using the kernels CS and M5. The number of neighbors is set to 100 for the M5 runs and 50 for the CS ones. To assess the convergence properties of the radial profiles of the considered hydrodynamic variables, the M5 run was replicated using 200   000 particles (M5 − 200 k). All of these runs were performed by incorporating in the SPH equations the AC term (18). For the CS kernel, a reference run with the AC term disabled was considered (CS-NOAC). The gravitational softening length was assumed to be εg = 0.02.

For each test case, we performed simulations with different settings of the AV parameter ld. However, introducing a time-dependent AV scheme in SPH reduces the severity of the numerical viscosity effects and in particular, for the test investigated here, produces a radial entropy profile at the shock front in better agreement with the 1D PPM reference solution. For this reason, the results for different test cases are presented for runs where a standard AV formulation was used (runs AV0 of V11). This is done to highlight differences in the simulation results due to the use of different kernels and the inclusion of the AC term in the SPH thermal energy equation. This choice of simulation parameters allows one to compare the test calculations performed here with previous results presented in Sect. 5.2 of V11.

The average radial profiles at t = 0.8 of density, entropy, time-dependent AC parameter, and radial velocity are displayed in Fig. 12. The black solid lines indicate the profiles of the 1D PPM calculation of Steinmetz & Mueller (1993). Broadly speaking, one expects to see some differences between the profiles of the NOAC simulation and those of the other runs. This should be valid in particular for the entropy profiles. However, at the considered epoch, Fig. 12 shows that the entropy profile of the NOAC run is still very similar to the others. This occurs because at t = 0.8 the shock front propagation is still in the early phase and the smoothing in the thermal energy due to the AC term is not yet significant.

thumbnail Fig. 13

As in Fig. 12, a closer view at t = 0.9 of the radial profiles of density and entropy in the proximity of the shock front.

Open with DEXTER

The behavior of the αC profiles is in accordance with that expected. In particular, the profiles exhibit a peak that corresponds to the shock front location, which occurs approximately at r ~ 0.18, and a decay as one moves away from it. The profile shapes depend on the kernel shape and weakly on resolution. The former is interpreted as a consequence of the improved mixing capability provided by the higher order of the kernel. We note that there is an increase in the profiles as one moves toward the sphere center, which occurs because of the presence of a temperature gradient. However, in Eq. (18) this radial dependence is of no consequence because behind the shock front the sphere quickly achieves a hydrostatic equilibrium and for r ≲ 0.1.

The epoch examined in Fig. 12 was chosen so as to compare the simulation profiles with the corresponding ones presented in Fig. 40 of Springel (2010b) and realized using the new moving-mesh code AREPO. In particular, a visual inspection shows that ahead of the shock front the accuracy of the entropy profiles shown in Fig. 12 can be considered to be in-between that of the profile produced using the AREPO code and the one realized with the SPH code Gadget-2 (right panels in Fig. 40 of Springel 2010b, test cases “moving-mesh” and “SPH”, respectively). We note, however, that a strict comparison between the different simulation profiles is difficult because the number of resolution elements employed by Springel (2010b) in the cold-gas sphere tests are lower by about  ~1/3 than the one used here. In particular, the smaller shock broadening found here can be ascribed to a resolution effect.

Finally, Fig. 13 shows a closer view of the radial density and entropy profiles in the proximity of the shock front at t = 0.9. The NOAC simulation exhibits the worst accord with the PPM reference solution, while for the high-resolution run M5 − 200k the improvement in accuracy is minimal with respect to the other AC runs. This suggests that to obtain a significant improvement in the adherence of the simulation profiles with the PPM solution one must use a much higher number of simulation elements (say  ≳ 106, Springel 2010b).

A noteworthy feature of the profiles of Fig. 13 is that now the AC simulation profiles are in much better agreement with the PPM profiles. This is interpreted as a numerical effect in which the pre-shock entropy, which is generated by the AV implementation ahead of the shock front owing to the converging flow, is now strongly reduced because the AC term in the SPH energy equation now removes this excess of internal energy. This in turn implies that, behind the shock front, a thermal energy profile is closer to the PPM solution profile.

Comparisons of the profiles of Fig. 13 with the corresponding ones displayed in Fig. 4 of V11 show that the dependence of the simulation profiles on the adopted AC scheme is stronger than that due to the time-dependent AV formulation. In fact, these results are almost unchanged when one performs SPH simulations using a time-dependent AV scheme in place of the standard one.

Finally, two runs were performed using the CRT and LIQ kernels and a number of neighbors set to Nsph = 50. The profiles of the two simulations have not been shown here to avoid overcrowding in the panels. The performances of the CRT run are quite similar to those of the CS one, whilst the profiles of the LIQ simulation are largely inconsistent with the reference PPM solution profiles.

4.2. 2D Rayleigh-Taylor instability

The Rayleigh-Taylor (RT) instability arises when a heavier fluid is placed on top of a lighter fluid (Chandrasekhar 1961) in the presence of an external gravitational field. The fluids are in pressure equilibrium with the external field, and in this configuration the system is unstable in the presence of small perturbations at the interface. The lighter fluid will then begin to rise and the denser one to fall. This process leads to the development of characteristic finger-like structures before the fluids enter the non-linear phase where they completely mix together. To evaluate the ability of the AC-SPH code to properly describe the evolution of RT instabilities, we consider a 2D version of the test. The initial conditions were chosen to be similar to the numerical test implemented by Abell (2011) to validate its new rpSPH formulation of SPH equations, so as to compare Abell’s results consistently. The computational domain consists in a 2D box with coordinates x ∈  {0,1/2}, y ∈  { 0,1 } . The boundaries are periodic in x and reflecting in y. The density is ρ1 = 2 at the top and ρ2 = 1 at the bottom, with a density profile (48)where Δy = 0.05. This ensures a smoothing in density at the interface y = 0.5 that allows a consistent numerical behavior (Robertson et al. 2010). This density profile is realized by constructing a HCP lattice of N = 6202 equal mass particles in which the spacing is varied according to the procedures described in Sect. 3.1.1 until the profile of Eq. (48) is satisfied.

thumbnail Fig. 14

Density maps of the 2D Rayleigh-Taylor tests are shown at time t = 3 for standard SPH (NOAC) and AC-SPH simulations performed using N = 6202 particles with different kernels and neighbor numbers. The density distributions can be compared directly with the map shown in Fig. 4 of (Abell 2011,top right panel), as the adopted initial conditions for the tests are the same.

Open with DEXTER

The pressure at the interface is set to P0 = ρ1/γ = 10/7, where γ = 1.4 and varies with y according to P(y) = P0 − (y)(y − 1/2), g = 1/2 being the external acceleration, so that the system is initially in hydrostatic equilibrium. For the particles that satisfy the condition  | y − 0.5 |  < 0.2, a velocity perturbation is applied in the y direction given by (49)where δvy = 0.1.

We present results for the CS, CRT, M5, and, for reasons that we discuss later, we consider also the M6 kernel (Price 2012a, Sect. 2.3).

The number of neighbors for these runs is set to Nsph = 32, 32, 50, 110, respectively. As in the other test cases considered here, for comparative purposes, we ran a simulation (CS-NOAC) in which the AC term in the thermal energy equation was disabled. The following settings were used for the AV parameters  { αmin,αmax,ld }  =  { 0.1,1.5,0.2 } .

In Fig. 14, at the time t = 0.3, we present the 2D density maps of the different RT tests. As expected, standard SPH is unable to correctly capture the development of the RT instability, as indicated by the leftmost panel of Fig. 14 (CS kernel, NOAC run). However, from this viewpoint the improvement is minimal even for the corresponding AC version of the considered run (CS-32). On the contrary, a significant improvement is obtained if one uses the M5 kernel or the next in order M6, which is a quintic polynomial (Price 2012a). In fact, in the rightmost panel of Fig. 14 (M6 run), the center of mass of the RT spikes appears at a lower position than in the M5 run and therefore the convergence might still not have been reached yet, even for the M6 run.

This strong dependence of the code performances on the kernel order demonstrates that, for the RT test, accuracy in pressure forces is a fundamental issue. These results are consistent with those of Sect. 3.1 and indicate that the poor performances of SPH when handling KH or RT instabilities, in particular for very subsonic flows, are mainly due to the leading errors in the momentum equation. These errors decrease as the order of the kernel increases, which implies that the accuracy in pressure force estimates is higher and thus the velocity noise is lower. This, in turn, indicates a better capability to capture the growth of the instability.

Similar results were obtained by McNally et al. (2012), who ran a suite of 2D KH simulations using carefully crafted initial conditions with the goal of assessing different hydrodynamic code capabilities. The authors conclude that for SPH the code performances are strongly related to the order of the kernel employed in the simulation and thus to the accuracy of the gradient estimates.

Moreover, the simulation behavior of Fig. 14 is consistent with the findings of García-Senz et al. (2012), who present a new formulation of SPH in which gradients are estimated using a tensorial approach that conserves both linear and angular momentum. The authors demonstrate that when using this higher-order gradient estimator the interpolation of physical quantities is significantly more accurate than in the standard method. In particular, several tests showed that the code is now able to successfully capture the development of KH and RT instabilities.

To summarize, the results of this section demonstrate that the performance of SPH in modeling hydrodynamic subsonic instabilities depend critically on the accuracy of the gradient estimator. The formulation proposed by García-Senz et al. (2012) looks particularly promising in this aspect in overcoming the present SPH difficulties.

Finally, it must be stressed that the simulation results of Fig. 14 were obtained by setting δvy = 0.1. If the same test had been performed using δvy = 0.01, as in the RT test of the new moving-mesh code AREPO (Springel 2010b, Sect. 8.8), the velocity noise would have suppressed the growth of the instability even for the M6 run. We note that in their RT test, García-Senz et al. (2012) were able to recover the growth of the instability using the same amplitude of the velocity perturbation, δvy = 0.01.

4.3. Cluster comparison

We now discuss a set of cosmological cluster simulations that we constructed using the standard SPH formulation as well as the new AC-SPH version. We only considered non-radiative, or “adiabatic” simulations, in which the hydrodynamics were modeled according to the formulation presented in Sect. 2. Previous investigations (Frenk et al. 1999; Wadsley et al. 2004, 2008; Mitchell et al. 2009; Springel 2010b) showed that for the same initial conditions there are systematic differences between results extracted from cluster simulations produced using standard SPH and Eulerian AMR codes. Specifically, the level of central entropy is found to be lower in SPH simulations than in the corresponding AMR runs. The latter simulations are characterized by the presence of a flat entropy core, whereas the radial entropy profile produced by SPH runs increases steadily with radius. As outlined in the Introduction, the origin of this discrepancy has been the subject of intense debate (Agertz et al. 2007; Mitchell et al. 2009; Read et al. 2010) that has given as the main explanation the different degrees of numerical mixing present in the two codes.

To assess the capability of the AC-SPH code in solving this issue, we then ran several cluster simulations. A comprehensive study of the differences between the hydrodynamic variables of cluster simulated using the standard and the AC version of the SPH code will be presented in a forthcoming paper. Here, we only show the final radial entropy profile of the simulated clusters as the main variable that can be used to test the effectiveness of the AC-SPH formulation in solving the entropy problem in cluster cores.

A detailed description of the procedures used to construct the simulated cluster sample is given in Sect. 2 of V11, and we provide a brief summary here. The simulations were carried out assuming a spatially flat ΛCDM model, with matter density parameter Ωm = 0.3, vacuum energy density ΩΛ = 0.7, baryonic density Ωb = 0.0486, and h = H/100   km   s-1   Mpc-1 = 0.7. The scale-invariant power spectrum is normalized to σ8 = 0.9 on an 8   h-1 Mpc scale at the present age t0. An N-body cosmological simulation involving only gravity was first run with a comoving box of size L2 = 200   h-1 Mpc, starting from an initial time tin and ending at the final epoch t0. At this epoch, clusters of galaxies were identified in the simulation box as groups of particles that are associated with overdensities approximately in excess of . Several of these clusters were then resimulated individually using the AC-SPH code described here, coupled with a treecode gravity solver.

Initial conditions for a specified cluster were generated as follows: a spherical region with origin located at the cluster center was populated with a high-resolution (HR) grid, the radius of the HR sphere being such that it contains all of the original particles identified at t = t0 as cluster members. Both a gas particle and a dark matter particle were associated with each grid node, whose positions were perturbed according to the random realization of the original cosmological simulations and only particles located within the HR sphere were kept for the hydro run. The HR particles were surrounded by a low-resolution shell of dark matter particles, extracted from a grid with spacing twice that of the HR grid, the shell being introduced for the purpose of mimicking the effects of tidal forces. The grid spacing was chosen such that at the end of the procedure a cluster was simulated with Ngas ~ 220   000 gas particles and Ndm ~ Ngas dark matter particles in the inner HR sphere, whereas were used in the low-resolution shell. Particle masses were assigned according to the values of Ωm and Ωb.

The clusters selected to be re-simulated hydrodynamically were chosen with the following criterion. Originally (V11), the procedure previously described was repeated two more times to construct two new cluster samples extracted from cosmological simulations with box sizes L4 = 2L2 and L8 = 2L4. The three cluster samples were then combined to construct a final sample Sall of  ~ 160 clusters covering nearly a decade in cluster masses. Cluster members of sample Sall were then ranked according to their dynamical state, as measured by an appropriate statistical indicator. Two sub-samples of four test clusters each, denoted respectively by Q and P, were then constructed by extracting from sample Sall clusters with membership criterion for sample Q (P) of being in a fully relaxed (perturbed) state. These two sub-samples were then used to investigate the effects of AV in hydrodynamic SPH simulations of galaxy clusters (V11). Here the four test clusters of sub-sample Q were chosen to investigate the performances of the AC-SPH code. This choice was motivated by the need to proper disentangle the amount of entropy mixing produced during merger processes from that specific to the numerical method in the final cluster entropy. The clusters of sub-sample Q were carefully selected on the basis of their dynamical state and are among the most relaxed of sample Sall. Therefore, for these clusters, differences in the final radial entropy profile between standard and AC runs can be safely ascribed to the new AC term implemented in the code. The AV simulation parameters of the simulations were  { αmin,αmax,ld }  =  { 0.1,1.5,0.2 } .

thumbnail Fig. 15

Final radial entropy profiles as a function of r/r200 for the four relaxed test clusters. The gas entropy is plotted in units of S200 Different line styles refer to different clusters, the numbers indicating the cluster membership in the cosmological ensemble from which they have been extracted. Thin (black) lines are for the profiles of runs performed using the standard SPH formulation, while thick (red) lines refer to the AC-SPH runs incorporating the AC term.

Open with DEXTER

The final entropy profiles for the four test clusters are displayed in Fig. 15, where the cluster entropies have been normalized to (50)where fb = Ωbm is the global baryon fraction and denotes the mass contained in a sphere of radius rΔ with mean density Δ times the critical density ρc = 3H2/8πG. As is customary, the cluster mass is defined by setting Δ = 200 and r200 denotes the cluster radius.

The results of Fig. 15 indicate that, although there is a certain degree of scatter between the entropy profiles of individual clusters, all of the AC-SPH runs consistently exhibit much shallower entropy profiles and higher core entropies than in the standard SPH runs at r/r200 ≲ 0.1. For a given cluster, the difference in the levels of central entropies produced by the two codes is about a factor  ~4, and now the values of central entropies are comparable with those produced using AMR codes in cosmological cluster simulations (Voit et al. 2005) and idealized cluster binary mergers (Mitchell et al. 2009). It should, however, be stressed that in Eulerian hydrodynamics it is the numerical scheme that forces the fluid to be mixed below the minimum cell size. This suggests that AMR simulations of galaxy clusters might overestimate the correct level of core entropy because of fluid mixing (Springel 2010b). To summarize, the development of entropy cores in the AC-SPH runs is clearly a consequence of the heat diffusion term described by Eq. (18) now present in the energy equation. This term acts to redistribute the internal energy produced in the shocks during dissipative processes, so that the subsequent entropy mixing leads, in the end, to higher levels of core entropies.

Similar results were obtained by Wadsley et al. (2008), who similarly added a heat diffusion term to the SPH energy equation, but using the prescription of Eq. (26) instead of the Eq. (18) adopted here. According to Wadsley et al. (2008), the value of the coefficient C tends to be problem dependent and for the cluster simulations the best agreement with the entropy profiles extracted from the mesh code runs is obtained by setting C ≃ 0.1−1. Here, the time-dependent formulation (22) presents the advantage of minimizing thermal diffusion away from shocks, as demonstrated in the cold gas sphere test by Fig. 12, where the radial profile of the AC parameter αAC is peaked at the shock location.

The simulations presented in this section were realized using adiabatic gas physics. However, a realistic modeling of the intracluster medium physics must incorporate the possibility that the gas cools radiatively. Previous simulations (V11) showed that the presence of cooling leads to the development in the inner cluster regions of dense compact cool-gas cores and subsequently of high levels of turbulence, the latter being produced by the hydrodynamical instabilities generated by the interaction of the compact cool cores with the ambient medium.

How this scenario is affected by incorporating into the SPH equations a numerical heat diffusion term is an issue that can be properly addressed only by resorting to numerical simulations. Here, we outline that two competing effects are expected to influence the final level of turbulence in a cluster core: the improved capability of the code to describe the development of hydrodynamical instabilities, and thus an increase in the amount of turbulence, and the reduced availability of cool gas because of the presence of material of higher entropy, which has been brought in the inner core by the enhanced fluid-mixing properties of the code.

Finally, a resolution study previously performed in V11 found for the simulation resolution employed here that there is a substantial stability in the entropy profiles of the simulations. Although the simulations of V11 do not incorporate the AC term, we do not expect strong variations in the resolution dependence of the profiles, given the similarity in the level of core entropies with those found using mesh-based codes.

5. Summary and conclusions

We have presented an SPH numerical scheme that incorporates an artificial conductivity term and uses an appropriate signal velocity for simulations including gravity. The AC formulation was introduced by Price (2008) as a solution to the problems encountered by standard SPH to correctly follow the development of KH instabilities, owing to the inconsistencies of the standard formulation in the description of density at a contact discontinuity. Here, a suite of hydrodynamic test problems has been investigated with the purpose of validating the new AC-SPH code and assessing its performances when using the specifically adopted signal velocity.

The results of the KH instability test were presented in Sect. 3.1, in which the code capabilities were tested by considering SPH simulations of the KH test performed using a large variety of initial condition set-up and SPH kernels. The set of initial conditions was chosen as in Valcke et al. (2010), so as to consistently compare the results. These are in accordance with the corresponding ones of Valcke et al. (2010) and indicate that, for the version of the KH test analyzed here, the AC implementation is important for the long-term behavior of the simulation.

Moreover, the results of the KH test indicate (Valcke et al. 2010; Read et al. 2010; McNally et al. 2012) that the poor performances of standard SPH in properly treating KH instabilities can be explained in terms of two distinct effects. The first is a general problem of consistency, which for the problem under consideration requires that smooth interfaces should be present at contact discontinuities, in order to obtain numerical convergence. The second effect, which in standard SPH suppresses the growth of KH instabilities, is the leading error in the momentum equation, which is caused by incomplete kernel sampling (Read et al. 2010) and quantified by the norm E0 defined by Eq. (36).

To circumvent this problem, several proposals have been made (Valcke et al. 2010; Read et al. 2010) in which the standard SPH cubic spline M4 kernel is replaced by the new LIQ or CRT kernels with steeper central profiles. These kernels present the advantage of being stable against particle clumping, so that the number of neighbors can be safely increased to reduce the error E0.

In this paper, we have considered the possibility of reducing sampling errors by considering higher-order B-spline kernels, specifically the M5 or quartic spline kernel. A striking result of the KH runs presented in Sect. 3.1 is that, for a given value of the ratio η of the smoothing length to the mean interparticle spacing, simulations of the KH test performed using the M5 kernel have amplitudes of the E0 error substantially smaller and in line with the behavior of the same quantity for simulations employing either the LIQ or CRT kernels. A linear stability analysis reveals that this result follows owing to the very good stability properties of the M5 kernel. The analysis indeed suggests that the clumping instability is absent for values of η up to η ≲ 2, which in 3D corresponds to Nsph ~ 520 neighbors.

These findings are consistent with the results of Dehnen & Aly (2012), who proposed to adopt, with the specific purpose of avoiding the clumping instability, the Wedland (1995) functions as a new class of kernels. The authors showed that in terms of stability and accuracy properties, the quartic spline performs extremely well when compared with the proposed Wedland functions. These results are strictly connected to the properties of the kernel Fourier transform, which according to the authors must be non-negative to avoid particle clumping, and are consistent with the findings of Sect. 3.1.3 since when increasing the kernel order both the B-spline and the Wedland functions approach the Gaussian.

We therefore propose, as a compromise between the need to reduce sampling errors while keeping the computational cost to a minimum, the use of the M5 kernel with a neighbor number in the range Nsph ~ 60−120 as the standard combination that guarantees sufficient accuracy in many SPH simulations of astrophysical problems. Moreover, the results of the Sod shock-tube test of Sect. 3.2 demonstrate that to obtain simulation profiles in accordance with the analytic solution, for simulations employing kernels with a modified shape the use of a much larger number of neighbors than in the case of the M5 runs is necessary.

The results of the gravity tests show that the adoption of the AC-SPH scheme significantly reduces, at the level tested in this paper, the differences seen in the hydrodynamics between standard SPH and grid-based simulations of self-gravitating structures. For the cold gas sphere, the entropy profile is in closer agreement with the PPM reference solution and for the cosmological cluster simulations, a key result is the final level of core entropies, which are consistent with those of the central entropies produced using AMR codes. Thus, it appears that in hydrodynamic simulations where self-gravity is important the AC term, accompanied by the proposed signal velocity, plays a key role as a mechanism for redistributing thermal energy and hence as a source of entropy mixing.

To summarize, results extracted from simulations of hydrodynamic tests where self-gravity is the dominant factor are in much better agreement with the corresponding ones obtained using mesh-based codes. These results then demonstrate the capability of the implemented AC-SPH scheme in properly following the formation of cosmic structures.

We note that the artificial heat conduction term was originally proposed by Price (2008) with the purpose of avoiding the inconsistencies encountered by standard SPH in the presence of density steps at contact discontinuities. A complementary view was proposed by Wadsley et al. (2008), who introduced the same term, albeit in a different numerical formulation, with the aim of modeling the level of diffusion due to turbulence. The two interpretations are not mutually inconsistent, although the results of the self-gravity tests presented in this paper support the view of a heat diffusion term that in SPH is capable of mimicking the diffusion due to turbulence. In a similar fashion, Violeau & Issa (2007) presented an SPH scheme to model a free-surface incompressible flow that, in analogy with 3D Large Eddy Simulations, assumes a Smagorinsky (1963) model for the filtered Navier-Stokes equations.

The results of the blob test simulations demonstrate that the instabilities leading to the expected cloud disruption can develop only when the SPH energy equation incorporates the AC term. A particularly interesting result is that an appropriate limiting condition must be implemented on the AC coefficients αC in order to avoid an unphysical amount of heat diffusion, which in turn leads to a cloud disruption that occurs too early. This limiter was identified with the Prandtl number and, for the AC signal velocity adopted here, severely limits the amplitude of the AC coefficients in the regime of strong supersonic flows.

The AC-SPH simulations of the blob test incorporating now the new constraint support this view, since Fig. 11 shows for the new runs a cloud mass-loss rate that is in closer agreement with the rates obtained from simulations realized using a completely independent numerical scheme (Springel 2011). However, it must be stressed that the condition specified by Eq. (43) was calculated for a perfect monoatomic gas with γ = 5/3, so that a physically motivated constraint on the artificial heat diffusion of the simulated medium should be considered problem-dependent.

However, the code still has several problems that render its use problematic if the development of hydrodynamic instabilities need to be followed in the regime of subsonic flows. The results of the simulations indicate that these shortcomings are not due to the AC implementation, but rather are intrinsic to the standard formulation with which gradients are calculated in SPH and the related errors are subsequently introduced in the momentum equation. In particular, simulations of the RT test show that increasing the kernel order alleviates the problems but does not solve them. These results are in accordance with recent findings (McNally et al. 2012; García-Senz et al. 2012) and clearly demonstrate that for very subsonic flows, the poor performances of SPH in modeling hydrodynamic instabilities are strictly connected to the code accuracy in gradient estimates.

The formulation of García-Senz et al. (2012), which was proposed with the aim of calculating SPH gradients with accuracy as high as possible while keeping the benefits of a Lagrangian formulation, looks in this respect very promising and encourages further investigations along the numerical approach proposed in this paper. These would be particularly relevant in the light of the results of Bauer & Springel (2012), who claim that standard SPH fails to properly model the regime of subsonic turbulence. They reach their conclusions by comparing results extracted from simulations of driven subsonic turbulence realized using the moving-mesh code AREPO and GADGET SPH. Their results were criticized by Price (2012a), for whom the use of a time-dependent AV is critical in SPH simulations of subsonic turbulence.

A re-analysis of the Bauer & Springel (2012) simulations using the AC-SPH code presented here, augmented with improved gradient operators, would then be fundamental for achieving a deeper understanding of the capability of the different numerical methods to model subsonic turbulence, which is expected to have a significant impact on shaping the thermodynamic properties of baryons in cosmological haloes and, subsequently, the process of galaxy formation.


1

Note that there is a misprint in the sign of the corresponding Eq. (28) of Price (2008).

References

  1. Abel, T. 2011, MNRAS, 413, 271 [NASA ADS] [CrossRef] [Google Scholar]
  2. Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963 [NASA ADS] [CrossRef] [Google Scholar]
  3. Balsara, D. 1995, J. Comp. Phys., 121, 357 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bauer, A., & Springel, V. 2012, MNRAS, 423, 2558 [NASA ADS] [CrossRef] [Google Scholar]
  5. Berger, M. J., & Colella, P. 1989, J. Comp. Phys., 82, 64 [NASA ADS] [CrossRef] [Google Scholar]
  6. Blundell, S. J., & Blundell, K. M. 2006, Concepts in Thermal Physics (Oxford University Press) [Google Scholar]
  7. Børve, S., Omang, M., & Trulsen, J. 2004, ApJS, 153, 447 [NASA ADS] [CrossRef] [Google Scholar]
  8. Brookshaw, L. 1985, Proc. Astron. Soc. Australia, 6, 207 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cha, S.-H., Inutsuka, S.-I., & Nayakshin, S. 2010, MNRAS, 403, 1165 [NASA ADS] [CrossRef] [Google Scholar]
  10. Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford University Press) [Google Scholar]
  11. Colella, P., & Woodward, P. R. 1984, J. Comp. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  12. Cullen, L., & Dehnen, W. 2011, MNRAS, 408, 669 [NASA ADS] [CrossRef] [Google Scholar]
  13. Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dimotakis, P. E. 2005, Ann. Rev. Fluid Mech., 37, 329 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  15. Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753 [NASA ADS] [CrossRef] [Google Scholar]
  16. Evrard, A. E. 1998, MNRAS, 235, 911 [NASA ADS] [CrossRef] [Google Scholar]
  17. Frenk, C. S., White, S. D. M., Bode, P., et al. 1999, ApJ, 525, 554 [NASA ADS] [CrossRef] [Google Scholar]
  18. Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
  19. Fulk, D. A., & Quinn, D. W. 1996, J. Comp. Phys., 126, 165 [NASA ADS] [CrossRef] [Google Scholar]
  20. García-Senz, D., Cabezón, R. M., & Escartín, J. A. 2012, A&A, 538, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375 [NASA ADS] [CrossRef] [Google Scholar]
  22. Hernquist, L., & Katz, N. 1989, ApJS, 70, 419 [NASA ADS] [CrossRef] [Google Scholar]
  23. Heß, S., & Springel, V. 2010, MNRAS, 406, 2289 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hubber, D. A., Batty, C. P., McLeod, A., & Whitworth, A. P. 2011, A&A, 529, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Inutsuka, S.-I. 2002, J. Comp. Phys., 179, 238 [NASA ADS] [CrossRef] [Google Scholar]
  26. Junk, V., Walch, S., Heitsch, F., et al. 2010, MNRAS, 407, 1933 [NASA ADS] [CrossRef] [Google Scholar]
  27. Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73 [NASA ADS] [CrossRef] [Google Scholar]
  28. Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212 [NASA ADS] [Google Scholar]
  29. Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] [Google Scholar]
  30. McNally, C. P., Lyra, W., & Passy, J.-C. 2012, ApJS, 201, 18 [NASA ADS] [CrossRef] [Google Scholar]
  31. Merlin, E., Buonomo, U., Grassi, T., Piovan, L., & Chiosi, C. 2010, A&A, 513, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Mitchell, N. L., McCarthy, I. G., Bower, R. G., Theuns, T., & Crain, R. A. 2009, MNRAS, 395, 180 [NASA ADS] [CrossRef] [Google Scholar]
  33. Monaghan, J. J. 1997, J. Comp. Phys., 136, 298 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  34. Morris, J. P., & Monaghan, J. J. 1997, J. Comp. Phys., 136, 41 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  35. Monaghan, J. J. 2000, J. Comp. Phys., 159, 290 [NASA ADS] [CrossRef] [Google Scholar]
  36. Monaghan, J. J. 2005, Rep. Progress Phys., 68, 1703 [NASA ADS] [CrossRef] [Google Scholar]
  37. Morris, J. P. 1996, Analysis of Smoothed Particle Hydrodynamics with Applications, Ph.D. Thesis, Monash University, Melbourne [Google Scholar]
  38. Murante, G., Borgani, S., Brunino, R., & Cha, S.-H. 2011, MNRAS, 417, 136 [NASA ADS] [CrossRef] [Google Scholar]
  39. Norman, M. L. 2005, The Impact of AMR in Numerical Astrophysics and Cosmology, in Adaptive Mesh Refinement – Theory and Applications, eds. T. Plewa, T. Linde, & V. G. Weirs (Berlin, New York: Springer), Lect. Notes Comput. Sci. Eng., 41, 413 [Google Scholar]
  40. Norman, M. L., & Bryan, G. L. 1999, in Numerical Astrophysics, eds. S. M. Miyama, K. Tomisaka, & T. Hanawa (Kluwer, Boston), Astrophys. Space Sci. Lib., 240, 19 [Google Scholar]
  41. O’Shea, B. W., Bryan, G., Bordner, J., et al. 2005, Introducing Enzo, an AMR Cosmology Application, in Adaptive Mesh Refinement – Theory and Applications, eds. T. Plewa, T. Linde, & V. G. Weirs (Berlin, New York: Springer), Lect. Notes Comput. Sci. Eng., 41, 341 [Google Scholar]
  42. Price, D. J. 2005, Ph.D. Thesis [arXiv:astro-ph/0507472] [Google Scholar]
  43. Price, D. J. 2008, J. Comp. Phys., 227, 10040 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  44. Price, D. J. 2012a, J. Comp. Phys., 231, 759 [NASA ADS] [CrossRef] [Google Scholar]
  45. Price, D. J. 2012b, MNRAS, 402, L33 [NASA ADS] [CrossRef] [Google Scholar]
  46. Price, D. J., & Monaghan, J. J. 2005, MNRAS, 364, 384 [NASA ADS] [Google Scholar]
  47. Rasio, F. A., & Lombardi, J. C. Jr. 1999, J. Comp. Phys., 109, 213 [Google Scholar]
  48. Robertson, B. E., Kravtsov, A. V., Gnedin, N. Y., Abel, T., & Rudd, D. H. 2010, MNRAS, 401, 2463 [NASA ADS] [CrossRef] [Google Scholar]
  49. Read, J. I., & Hayfield, T. 2011, MNRAS, 422, 3037 [NASA ADS] [CrossRef] [Google Scholar]
  50. Read, J. I., Hayfield, T., & Agertz, O. 2010, MNRAS, 405, 1513 [NASA ADS] [Google Scholar]
  51. Rosswog, S. 2009, New Astron. Rev., 53, 78 [NASA ADS] [CrossRef] [Google Scholar]
  52. Rosswog, S., & Price, D. J. 2007, MNRAS, 379, 915 [NASA ADS] [CrossRef] [Google Scholar]
  53. Rosswog, S., Davies, M. B., Thielemann, F.-K., & Piran, T. 2000, A&A, 360, 171 [NASA ADS] [Google Scholar]
  54. Ritchie, B. W., & Thomas, P. A. 1991, MNRAS, 323, 743 [NASA ADS] [CrossRef] [Google Scholar]
  55. Ryu, D., Ostriker, J. P., Kang, H., & Cen, R. 1993, ApJ, 414, 1 [NASA ADS] [CrossRef] [Google Scholar]
  56. Saitoh, T. R., & Makino, J. 2009, ApJ, 697, L99 [NASA ADS] [CrossRef] [Google Scholar]
  57. Saitoh, T. R., & Makino, J. 2012, ApJ, submitted [arXiv:1202.4277] [Google Scholar]
  58. Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic Press) [Google Scholar]
  59. Sijacki, D., Vogelsberger, M., Keres, D., Springel, V., & Hernquist, L. 2012, MNRAS, 424, 2999 [NASA ADS] [CrossRef] [Google Scholar]
  60. Smagorinsky, J. 1963, Mon. Weather Rev., 91, 99 [NASA ADS] [CrossRef] [Google Scholar]
  61. Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
  62. Springel, V. 2010a, ARA&A, 48, 391 [NASA ADS] [CrossRef] [Google Scholar]
  63. Springel, V. 2010b, MNRAS, 401, 791 [NASA ADS] [CrossRef] [Google Scholar]
  64. Springel, V. 2011 [arXiv:1109.2218] [Google Scholar]
  65. Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649 [NASA ADS] [CrossRef] [Google Scholar]
  66. Steinmetz, M., & Mueller, E. 1993, A&A, 268, 391 [NASA ADS] [Google Scholar]
  67. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [NASA ADS] [CrossRef] [Google Scholar]
  68. Swegle, J. W., Hicks, D. L., & Attaway, S. W. 1995, J. Comp. Phys., 116, 123 [NASA ADS] [CrossRef] [Google Scholar]
  69. Tasker, E. J., Brunino, R., Mitchell, N. L., et al. 2008, MNRAS, 390, 1267 [NASA ADS] [CrossRef] [Google Scholar]
  70. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Valdarnini, R. 2011, A&A, 526, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Valcke, S., de Rijcke, S., Rödiger, E., & Dejonghe, H. 2010, MNRAS, 408, 71 [NASA ADS] [CrossRef] [Google Scholar]
  73. Violeau, D., & Issa, R. 2007, Int. J. Num. Meth. Fluids, 53, 277 [CrossRef] [Google Scholar]
  74. Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MNRAS, 364, 909 [NASA ADS] [CrossRef] [Google Scholar]
  75. Vogelsberger, M., Sijacki, D., Keres, D., Springel, V., & Hernquist, L. 2011 [arXiv:1109.1281] [Google Scholar]
  76. Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New Astron., 9, 137 [NASA ADS] [CrossRef] [Google Scholar]
  77. Wadsley, J. W., Veeravalli, G., & Couchman, H. M. P. 2008, MNRAS, 387, 427 [NASA ADS] [CrossRef] [Google Scholar]
  78. Wedland, H. 1995, Adv. Comp. Math., 4, 389 [CrossRef] [Google Scholar]
  79. Wetzstein, M., Nelson, A. F., Naab, T., & Burkert, A. 2009, ApJS, 184, 298 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

KH parameters for the simulations.

Table 2

Main simulation parameters of the KH tests.

Table 3

Average densities and sample standard deviations estimated from the SPH densities of a configuration of one million particles.

All Figures

thumbnail Fig. 1

Density maps for some of the 2D KH instability tests described in Sect. 3.1. From left to right, each column shows the panels for runs having initial conditions with the same Mach number at the corresponding time t = τKH, as listed in Table 1. From top to bottom, different rows refer to the simulations SPH, RHO, LIQ, and LP. The last one uses the linear quartic kernel but with the signal velocity (24), whereas the first three use the expression (25) (see Table 2). The plots can be compared directly with Fig. 10 of Valcke et al. (2010).

Open with DEXTER
In the text
thumbnail Fig. 2

As in Fig. 1 but for the set of simulations LIQ, CRT, and M5.

Open with DEXTER
In the text
thumbnail Fig. 3

Rendered plots of the αC parameters are shown for the same runs as in Fig. 2. The distribution of particle values  has been interpolated at the map grid points according to the SPH prescription.

Open with DEXTER
In the text
thumbnail Fig. 4

As in Figs. 1 and 2 but here at the time t = 2 for all of the panels. Note from Table 1 that this implies t ≫ τKH for runs with high Mach numbers.

Open with DEXTER
In the text
thumbnail Fig. 5

Time evolution of the velocity field amplitude in the y direction, as measured by the  = 12π mode of the Fourier transform of vy, for some of the KH instability tests described in Sect. 3.1. Each panel refers to KH simulations performed with the same Mach number, the initial conditions set-up being given in Table 1. Within each panel, the different curves are for AC-SPH runs with different simulation parameters, as specified in Table 2. The black solid line indicates the expected linear-theory growth rate.

Open with DEXTER
In the text
thumbnail Fig. 6

Top panels: distribution at time t = τKH of the errors  | Ei0 |  plotted versus y, as defined by Eq. (36), for the KH runs with Mach number M = 1 and different simulation parameters (See Table 2). The red histograms show the mean values of the binned distributions. Bottom panels: each panel shows the nearest neighbor map of the run in the corresponding top panel. This is defined as the distribution, interpolated at the map grid points of the normalized distances , where  is the distance  | xi − xk |  of the kth neighbor of the particle i and the neighbors are sorted so that .

Open with DEXTER
In the text
thumbnail Fig. 7

Contour plots of the dispersion relation given in Eq. (40) are shown as a function of wavenumber k and smoothing length h for a regular lattice of particles with spacing Δ, where we consider a wavevector with orientation k = k(1,0,0). Gray areas indicate the instability regions for which ω2 < 0. From the left to right, the top panels show in the case of the M5 kernel the instability regions of the longitudinal and the two transverse waves. The bottom panels refer to the CS (M4) kernel.

Open with DEXTER
In the text
thumbnail Fig. 8

Results at t = 0.12 of the 3D shock tube test. The profiles of density, pressure, thermal energy, and velocity are plotted clockwise from top left. The solid black line represents the analytical solution, while lines with different styles and colors are the profiles of the AC-SPH runs with different kernels and neighbor numbers, as illustrated in the pressure panel.

Open with DEXTER
In the text
thumbnail Fig. 9

Radial density profiles of the 3D Sedov blast wave test at t = 0.06. The solid black line indicates the analytic solution, while the simulation profiles are obtained by averaging for each radial bin SPH densities calculated from the particle distributions over a set of grid points located at the surface of a spherical shell and uniformly spaced in angular coordinates. All of the runs were performed using N = 524   288 particles.

Open with DEXTER
In the text
thumbnail Fig. 10

Density maps of the blob test in the central plane x = Lx/2 at t = 1, 2, 3 for SPH runs with (CS, M5) and without (CS NOAC) the AC term. Time is in units of τKH ~ 1.6τcr and the number accompanying the kernel label indicates the number of neighbors of the run. Axis units are in Mpc.

Open with DEXTER
In the text
thumbnail Fig. 11

Mass loss of the cloud versus time for the blob test simulations. Particles are defined as a cloud member if their temperatures and densities fulfill the conditions Ti < 0.9Text and ρi > 0.64ρcl. Left panel: simulations have been performed incorporating the AC term into the SPH equations and using the kernels CS, CRT, and M5. For the CS kernel, a simulation is also shown as reference without AC (NOAC, solid black line). Right panel: for the CS kernel, simulations with additional constraints on the particle AC parameters  have been performed using N = 106 (N6) and N = 105 (N5) particles, see text for more details. Black squares have been extracted from Fig. 10 of Springel (2011) and indicate the behavior of the cloud mass versus time in a similar test performed using the moving-mesh code AREPO and with (64 × 64 × 128) resolution elements.

Open with DEXTER
In the text
thumbnail Fig. 12

Radially averaged profiles at t = 0.8 of the Evrard collapse test. Clockwise from the top left panel: profiles of density, entropy, artificial conductivity parameter αC, and radial velocity. Curves with different line styles and colors refer to SPH runs performed using different kernels and with (AC) or without (NOAC) the AC term in the SPH energy equation. The black solid lines indicate the 1D PPM reference solution of Steinmetz & Mueller (1993).

Open with DEXTER
In the text
thumbnail Fig. 13

As in Fig. 12, a closer view at t = 0.9 of the radial profiles of density and entropy in the proximity of the shock front.

Open with DEXTER
In the text
thumbnail Fig. 14

Density maps of the 2D Rayleigh-Taylor tests are shown at time t = 3 for standard SPH (NOAC) and AC-SPH simulations performed using N = 6202 particles with different kernels and neighbor numbers. The density distributions can be compared directly with the map shown in Fig. 4 of (Abell 2011,top right panel), as the adopted initial conditions for the tests are the same.

Open with DEXTER
In the text
thumbnail Fig. 15

Final radial entropy profiles as a function of r/r200 for the four relaxed test clusters. The gas entropy is plotted in units of S200 Different line styles refer to different clusters, the numbers indicating the cluster membership in the cosmological ensemble from which they have been extracted. Thin (black) lines are for the profiles of runs performed using the standard SPH formulation, while thick (red) lines refer to the AC-SPH runs incorporating the AC term.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.