Gasoline2: a modern smoothed particle hydrodynamics code

Abstract

The methods in the Gasoline2 smoothed particle hydrodynamics (SPH) code are described and tested. Gasoline2 is the most recent version of the Gasoline code for parallel hydrodynamics and gravity with identical hydrodynamics to the Changa code. As with other Modern SPH codes, we prevent sharp jumps in time-steps, use upgraded kernels and larger neighbour numbers and employ local viscosity limiters. Unique features in Gasoline2 include its Geometric Density Average Force expression, explicit Turbulent Diffusion terms and Gradient-Based shock detection to limit artificial viscosity. This last feature allows Gasoline2 to completely avoid artificial viscosity in non-shocking compressive flows. We present a suite of tests demonstrating the value of these features with the same code configuration and parameter choices used for production simulations.

hydrodynamics, methods: numerical

1 INTRODUCTION

We present the Gasoline2¹ smoothed particle hydrodynamics (SPH) code, which is the culmination of improvements over several years to the original Gasoline SPH code (Wadsley, Stadel & Quinn 2004). The purpose of this paper is to provide details of the modern SPH implementation that is current in both Gasoline2 and the Changa code (Menon et al. 2015). Gasoline2 includes several unique features that substantially enhance the accuracy of SPH on problems of astrophysical interest. We explicitly note that several of these features were introduced previously (e.g. Wadsley, Veeravalli & Couchman 2008; Shen, Wadsley & Stinson 2010; Keller et al. 2014). The current work is the first to provide complete details and show results on standard tests when all these features are combined.

SPH is a method for discretizing and solving Euler's equations; introduced by two independent groups in Lucy (1977) and Gingold & Monaghan (1977). SPH discretizes comoving mass elements, making it a Lagrangian method. This has the advantage of making the tracing of fluid elements trivial. This also makes SPH automatically concentrate resolution elements where densities are largest, which is a common strategy for adaptive codes in astrophysics. SPH has perfect Galilean invariance and exact momentum, angular momentum and energy conservation due to its pairwise interactions. In order to handle shocks, SPH requires the use of artificial viscosity, typically paired with a limiter (e.g. Balsara 1995) to reduce the viscosity away from shocks. This makes SPH quite numerically robust compared to Riemann-solver based codes. SPH's particle nature is also ideal for unification with N-body schemes for gravity (Evrard 1988; Hernquist & Katz 1989). Lagrangian schemes benefit from time-steps that depend only on local wave speeds (e.g. sound speed) rather than absolute velocities making them efficient in cold, self-gravitating systems (e.g. discs) as well as avoiding frame and direction dependent behaviour.

SPH has become a pillar of numerical astrophysics, along with Eulerian adaptive-mesh refinement (AMR) schemes based on Riemann solvers such as enzo (Bryan & Norman 1997), flash (Fryxell et al. 2000), and ramses (Teyssier 2002). Recently these methods have been joined by hybrid schemes that attempt to merge the best features of both classes of method, including arepo (Springel 2010b) and gizmo (Hopkins 2015) (see also: Gaburov & Nitadori 2011).

While traditional SPH [as presented in Monaghan (1992, 2005), Gadget2 (Springel & Hernquist 2002), and the original Gasoline] has many attractive features, it became clear that it also had some significant shortcomings, especially when dealing with multiphase fluids. These problems were highlighted in Agertz et al. (2007) through the application of several major codes to a suit of novel test problems. Referring to SPH formalisms prior to ∼2008 as ‘traditional’ is fairly common in the literature (e.g. Gibson et al. 2009; Hopkins 2013). Other common labels include ‘standard’ SPH (Springel 2010a) and ‘classic’ SPH (Read, Hayfield & Agertz 2010; Bauer & Springel 2012). For the purposes of this work, we shall refer to all SPH in common usage prior to 2008 as ‘traditional’ SPH. We chose 2008 because it was the year that artificially suppressed mixing was first identified and then resolved for SPH (Wadsley et al. 2008; Price 2008), primarily in response to Agertz et al. (2007). Some of the issues raised were already known and SPH variants were proposed to address them (e.g. Ritchie & Thomas 2001), but were not widely adopted. It should also be noted that Gadget2 (Springel & Hernquist 2002) differed from other contemporary SPH codes in employing a formulation derived from a Lagrangian (first proposed by Nelson & Papaloizou 1994). This requires iterations to satisfy a constraint but offers specific benefits in simplifying entropy and energy conservation (as discussed further in Section 2). This change does not help resolve its primary shortcomings such as the mixing issue. In what follows, we present our own take on the state of SPH and important developments, focusing on those that influenced the design of Gasoline2. We refer the reader to recent reviews of SPH for additional perspectives, such as Springel (2010a); Price (2012a); and Rosswog (2015).

SPH and Eulerian codes generally agree well in simulations that involve supersonic flow (e.g. comparison projects by Kitsionas et al. 2009; Price & Federrath 2010; Bauer & Springel 2012), high-Mach number shocks (Tasker et al. 2008; Sijacki et al. 2012) and gravitational structure formation (O'Shea et al. 2005). However, in problems involving multiphase gas in near pressure equilibrium, SPH has been found to have a number of problems (e.g. Agertz et al. 2007; Read et al. 2010). These problems manifest as suppressed growth for fluid instabilities in shearing flows. Work by many authors, as described below, has demonstrated that these problems primarily arise from three sources: insufficient mixing, surface tension at density discontinuities and slow numerical convergence. The last issue is not about fundamental correctness but rather the amount of numerical effort required to get a satisfactory result.

A related issue is that formulations of SPH with constant artificial viscosity include a substantial shear viscosity term, which leads to excess angular momentum transport in differential rotation and subsonic turbulence (Bauer & Springel 2012). A popular partial solution has been the Balsara (1995) switch, but it is too noisy to strongly limit viscosity in general flows. An effective approach is time-dependent viscosity reduction away from shocks as first proposed by Morris & Monaghan (1997). We note that in the past many codes did not include this as a default option (e.g. Gadget2 and the original Gasoline). Bauer & Springel (2012) ran Gadget2 with a switch like that of Morris & Monaghan (1997) and achieved similar inertial ranges to grid codes (albeit with less small-scale structure). Thus, we do not include excess artificial viscosity as an outstanding problem for traditional SPH. Cullen & Dehnen (2010) (hereafter CD) proposed a more complex switch, which demonstrates little numerical angular momentum transport in various rotating tests. Price (2012b) and Hopkins (2013) have demonstrated excellent performance by modern SPH with time-dependent viscosity switches on shearing problems including subsonic turbulence.

The Lagrangian nature of SPH particles means that fluid quantities do not mix at small (sub-resolution) scales. This would be correct for laminar flows with limited molecular diffusion but most astrophysical flows have high Reynolds numbers and are expected to have a turbulent cascade that continues down below resolved scales. For example, Wadsley et al. (2008) showed that a lack of turbulent diffusion led to the low-entropy cores for SPH simulations of galaxy clusters in Frenk et al. (1999). Shen et al. (2010) showed that turbulent diffusion is required for realistic distributions of metals in and around simulated galaxies. These issues are most apparent in the presence of gravity, where buoyancy and convection separate particles with different entropies. Mesh codes diffuse all quantities due to numerical diffusion errors associated with advection, which can be close to the expected physical diffusion when the absolute fluid velocities are small. Problems arise in structured, high-Mach number flows, such as in typical disc galaxies (Benincasa et al. 2013).

A complete lack of mixing can suppress the growth of fluid instabilities, making traditional SPH unable to capture the growth of Kelvin–Helmholtz or Rayleigh–Taylor instabilities across discontinuous jumps in thermal energy (Price 2008; Read et al. 2010). This can lead to dense structures failing to mix within the correct time-scales (Agertz et al. 2007; Wadsley et al. 2008).

Every simulation of turbulent flows should include explicit turbulent transport terms, of which the leading order term is turbulent diffusion.

McNally, Lyra & Passy (2012) showed that explicit diffusion can be required to prevent numerical noise overtaking explicit perturbations such as in the Kelvin–Helmholtz instability. The most complex schemes explicitly evolve unresolved turbulent energy (e.g. Large-Eddy Schemes, Germano et al. 1991; Schmidt, Niemeyer & Hillebrandt 2006). For these reasons, in Gasoline2 we favour the idea of turbulent diffusion based on local velocity shear to mix fluid quantities. We note that approaches to mixing such as artificial conduction are also effective but can have side effects in gravitational fields, making an approach based on local velocity gradients more generally applicable (Price 2008).

Numerical surface tension exacerbates problems with mixing, realistic treatment of multiphase flow and two-phase instabilities. Surface tension arises from the formulation of pressure forces within traditional SPH and is related to inconsistencies in the pressure and force profiles at constant-pressure density jumps. Proposed solutions include estimating pressure directly rather than density (Ritchie & Thomas 2001) and higher numbers of neighbours (Read et al. 2010). The large number of neighbours and associated numerical expense prevented widespread adoption of these methods. Strong artificial conduction can reduce the effect of surface tension by removing sharp density gradients (Price 2008). Ritchie & Thomas (2001) promoted the idea of smoothing to calculate pressure directly and then deriving density from pressure. Hopkins (2013) employed pressure-entropy style SPH, combining this idea with an entropy evolution equation. He demonstrated that substantial improvement could be achieved with modified pressure-based formulations without large neighbour numbers. Gasoline2 uses a new pressure force formulation to achieve similar results.

In some regimes, SPH errors scale sub-linearly with resolution due to the ‘E0-error’, which is associated with the re-gridding noise that keeps SPH particles in a glass-like configuration in evolved flows (Read et al. 2010). The arepo method (Springel 2010b) helps resolve this by separating the motion of volume elements from that of the fluid at the cost of not being a strictly Lagrangian method. SPH can reduce this error with large neighbour numbers. However, traditional SPH kernel functions are unable to use neighbour numbers larger than N_smooth ∼ 50 due to a pairing instability, which effectively limits the resolution to N_smooth < 50. This issue was fully resolved by Dehnen & Aly (2012), who showed the origin of the pairing and that Wendland (1995) smoothing kernels do not suffer from it.

Some issues have only come to the fore due to recent enhancements in resolution. For example, highly resolved Sedov blasts (Supernova explosions) with temperature jumps exceeding factors of 100 have only been a feature of galaxy simulations in recent years. This led to the realization that hydrodynamics codes must avoid rapid spatial variations in individual time-steps (Saitoh & Makino 2009).

It should be noted that specific, idealized tests can still be problem free for specific subsets of these modern SPH features. Our goal has been to develop a code where a single set of parameters and algorithms gives good results on all tests relevant to our application areas.

Modern SPH codes include features that alleviate the issues present in traditional SPH. As noted earlier, limiting numerical viscosity is important and the CD approach is commonly considered the current state of the art. The use of new kernels rather than the traditional SPH M4 spline kernel allows for more neighbours and higher accuracy (e.g. the Gasoline derived OSPH; Read et al. 2010). Many codes use newer forms of the SPH equations of motion (e.g. the Gadget-derived PSPH; Hopkins 2013). Modifying pressure forces and using explicit mixing terms (e.g. gcd+; Kawata et al. 2013) both reduce surface tension problems. Many codes include a combination of these features (e.g. the Gadget derived codes of Hu et al. 2014 and Beck et al. 2016).

The paper presents a complete description of all the methods now present in Gasoline2 making it a state of the art, modern SPH code in Section 2. We demonstrate the accuracy and performance of the resulting code in Section 3. Gasoline has been continually upgraded, including turbulent diffusion from 2008 and geometric-average-density pressure forces (see Section 2.3) in 2011. These were sufficient to address major issues with fluid instabilities (as shown in Section 3). The primary feature that was lacking until recently was CD-like time-dependent viscosity. Gasoline2 also has unique features to limit viscosity in convergent flows, which previously published viscosity limiters fail to do. We provide detailed motivations for our preferred formulation and demonstrate the benefits with specific tests. We conclude by discussing future directions for Lagrangian hydrodynamics in Section 4.

1.1 Gasoline, Pkdgrav and Changa

All versions of the Gasoline code have been built upon the pkdgrav parallel tree code (Stadel 2001). Gasoline2 retains the core architecture of pkdgrav. This includes the KDK integration scheme, gravity solver and tree structure. In 2005, it was cutting edge with comparable or superior performance to other codes (e.g. Gadget2 and 3; Springel 2005). However, new programming approaches were needed to take advantage of systems with thousands of cores. The core gravity and parallel routines in Gasoline stopped active development around 2008 in response to the plans for new, scalable parallel alternatives. Pkdgrav continued its development as an N-body only code, culminating in the Pkdgrav3 code (Potter, Stadel & Teyssier 2016). Pkdgrav3 is optimized for solving gravity on large parallel systems with GPUs and has run over a trillion particles.

Combining both hydrodynamics (SPH) and gravity is challenging for parallelization with distributions of work and data highly clustered in both space (large density contrasts) and time (multiple time-steps). These distributions also differ for different species (e.g. dark matter, stars and gas). An effective approach is inherently parallel languages such as Charm++ (Kale & Krishnan 1993) based on c++. The Charm++ execution model incorporates functions as asynchronous, concurrent tasks, work migration and internal task scheduling. These are managed without the direct intervention of the programmer. This model allows simulations to adapt to work imbalance, node failures and heterogeneous architectures at run time (see also: Theuns et al. 2015). Changa is a TreeSPH code written in Charm++ that incorporates many of the gravity code improvements of pkdgrav and all of the SPH methods implemented in Gasoline2. We refer the reader to Menon et al. (2015) for details of code design and performance. The most important difference is the use of an oct-tree rather than a binary KD-tree for gravity. This changes the distribution of the gravity force errors and thus calculations involving gravity are not identical using Gasoline2 and Changa. The gravity improvements include fast multipole methodology (Potter et al. 2016), so that Changa is also intrinsically faster than Gasoline2 for N-body simulations. On the other hand, SPH is based on neighbour lists and it is thus possible to achieve identical SPH densities and forces to round off. Over time, even round off-level differences accumulate and lead to noticeable differences at the particle level.

Gasoline2 and Changa have been developed together with the objective of maintaining identical SPH methods and results. In most functions the code has been directly copied. We will focus on the current SPH methods employed in these codes in Section 2 and present results for standard test problems in Section 3. The results presented here were run with Gasoline2. Changa produces essentially identical results.

2 METHOD

2.1 SPH estimators, kernels and smoothing lengths

All SPH codes are based on local summations over neighbours. For some fluid quantity, f_j, known at particle positions, |$\boldsymbol {x}_j$|⁠, we get an SPH smoothed estimate as follows:

\begin{eqnarray} f_{i,{\rm smoothed}} = \sum _{j} f_j {\rm W}(\boldsymbol {x}_i-\boldsymbol {x}_j,h_i,h_j), \end{eqnarray}

(1)

where W is a general kernel function and h_j is a per-particle smoothing length, indicative of the range of interaction of particle j. It is standard practice to directly refer to quantities such as density that are derived this way with the smoothed qualifier. For momentum and energy conservation a symmetric expression is required in the force terms, which led many authors to symmetrize the kernel in all summations. Starting with Springel & Hernquist (2002) it became common to use un-symmetrized or gather estimates for non-force terms such as the density. For the density and other per-particle quantities, Gasoline2 uses the gather (i.e. one-sided) estimate, which is equivalent to

\begin{eqnarray} {\rm W}(\boldsymbol {x}_i-\boldsymbol {x}_j,h_i,h_j) = \frac{1}{h_i^3} W\left(\frac{r_{ij}}{h_i}\right), \end{eqnarray}

(2)

where |$\boldsymbol {r}_{ij} = \boldsymbol {x}_i-\boldsymbol {x}_j$| and r_ij is its magnitude. For the specific kernel function, W(q), we use the Wendland kernels Wendland (1995) that do not suffer from a pairing instability. We also follow Dehnen & Aly (2012) in adjusting the kernel weight at particle i to correct for the small bias in the density estimator. Gasoline2 can be run using other kernels such as the traditional cubic spline Monaghan (1992).

We define the smoothing length, h_i, as half the distance to the furthest neighbour, consistent with its use in the original SPH literature (e.g. Monaghan 1992). To find the smoothing length, we find an exact number of neighbours. Our standard neighbour number (as used on the tests below) is 200. Historically, many Gasoline simulations were run with 32–64 neighbours. Gasoline and Changa employ search algorithms based on priority queues that efficiently return complete one-sided neighbour lists and do not require iterations.

Gadget2 (Springel & Hernquist 2002) and many other codes weight neighbours in the same way as the density sum and require a small number of iterations to find h_i. Hopkins (2013) pointed out that the weightings used for neighbour selection constraints can be completely generalized independent of the density expression. A fixed neighbour count corresponds to a uniform weighting. It should be noted that some codes (e.g. Gadget2) do not allow large deviations from the typical number of neighbours. A fixed neighbour count constraint can be incorporated into SPH frameworks derived from a Lagrangian (Nelson & Papaloizou 1994) but such frameworks are not used here. Consequences and implications for the scheme are discussed in Section 2.3.

Density in Gasoline2 thus only depends on quantities known at particle i and the positions of its neighbours,

\begin{eqnarray} \rho _i &=\sum _{j} m_j W(r_{ij},h_i). \end{eqnarray}

(3)

This is computationally efficient, particularly when only a subset of particles requires updated densities.

2.2 SPH gradients

The one-sided gradient of the kernel is

\begin{eqnarray} \nabla _i {\rm W}(r_{ij},h_i) = (\boldsymbol {x}_i-\boldsymbol {x}_j) \frac{1}{h_i^5} W^{\prime }(q), \end{eqnarray}

(4)

where |$W^{\prime }(q) = \frac{1}{q} \frac{{\rm d}W}{{\rm d}q}$| and q = r_ij/h_i.

Most modern codes use local gradients in velocity to estimate diffusion and limit dissipation terms associated with artificial viscosity. These may be calculated at the same time as density via one-sided estimates. For example, the divergence and curl of the velocity field |$\boldsymbol {v}_i$| at the location of the particle can be constructed directly from components of the velocity gradient tensor,

\begin{eqnarray} {\boldsymbol V}_{\alpha \beta } = \frac{\mathrm{\partial} v_{\alpha }}{\mathrm{\partial} x_{\beta }}, \end{eqnarray}

(5)

where α and β run through the Cartesian axes.

We employ a straightforward estimator for the velocity gradient tensor is as follows:

\begin{eqnarray} {\boldsymbol V}_{\alpha \beta }\left. \right|_i &= \frac{\sum _{j} (v_{\alpha \,i}-v_{\alpha \,j}) (x_{\beta \,i}-x_{\beta \,j}) m_j \frac{1}{h_i^5} W^{\prime }\left(\frac{r_{ij}}{h_i}\right)}{\frac{1}{3} \sum _{j} r_{ij}^2 m_j \frac{1}{h_i^5} W^{\prime }\left(\frac{r_{ij}}{h_i}\right)}. \end{eqnarray}

(6)

The numerator is a standard SPH gradient estimate and the denominator expression provides the improvement. This estimator is exact in the case of uniform contraction or expansion and corresponds to CD equation (B1). CD point out that it is biased if there are significant density gradients. However, we find that this bias is modest compared to the effects of particle noise on the estimator. For the purposes of estimating artificial dissipation the simpler form is sufficient. In particular, we see no difference on test problems.

For symmetrized force terms, we require a symmetric gradient of the kernel, ∇_iW_ij,

\begin{eqnarray} \bar{\nabla _i W_{ij}}=\frac{1}{2} f_i \nabla _i W(r_{ij},h_i) + \frac{1}{2} f_j \nabla _j W(r_{ij},h_j). \end{eqnarray}

(7)

The term,

\begin{eqnarray} f_i= \frac{\sum _{j} \frac{m_j}{\rho _i} r_{ij}^2 W^{\prime }\left(\frac{r_{ij}}{h_i}\right)}{\sum _{j} \frac{m_j}{\rho _j} r_{ij}^2 W^{\prime }\left(\frac{r_{ij}}{h_i}\right)}, \end{eqnarray}

(8)

is a correction of order unity that can be compared to the f_i term in Springel & Hernquist (2002) and is discussed in Section 2.3.1. The symmetrized expression simply reverses sign on the interchange of i and j. In practice, we accumulate force terms via two separate summations to particles i and j with different smoothing lengths in each case.

2.3 Hydrodynamical equations

SPH is an approach to solving Euler's equations, given in Lagrangian form as

\begin{eqnarray} \frac{{\rm d}\,\rho }{{\rm d} t} = -\rho \nabla \cdot \boldsymbol {v}, \end{eqnarray}

(9)

\begin{eqnarray} \frac{{\rm d}\,\boldsymbol {v}}{{\rm d} t} = -\frac{\nabla P}{\rho } + \boldsymbol {g}, \end{eqnarray}

(10)

\begin{eqnarray} \frac{{\rm d}\,u}{{\rm d} t} = -\frac{P}{\rho } \nabla \cdot \boldsymbol {v} + \Gamma - \Lambda , \end{eqnarray}

(11)

where ρ, |$\boldsymbol {v}$| and u are the fluid density, velocity and internal energy, respectively, and the three equations express mass, momentum and energy conservation, respectively. |$\boldsymbol {g}$| represents body forces such as gravity. To close the system, we need an equation of state such as P = (γ − 1)ρu for an ideal gas, where γ is the ratio of the specific heats. The precise form depends on the complexity of the fluid being simulated and is closely connected to heating, Γ, and cooling, Λ, terms which are generally complex functions of density, temperature and composition. In what follows, we omit body forces and heating and cooling terms.

Equation (9) is automatically satisfied in SPH by equation (3).

The momentum equation is expressed as

\begin{eqnarray} \frac{{\rm d}\,\boldsymbol {v}_i}{{\rm d}t}& = & -\sum _{j} m_j \left({\frac{P_i+P_j}{\rho _i \rho _j}+\Pi _{ij}}\right) \bar{\nabla _i W_{ij}}, \end{eqnarray}

(12)

where P_j is pressure, |$\boldsymbol {v}_i$| velocity and Π_ij is an artificial viscosity term, discussed in Section 2.4. This Geometric Density Average Force (GDF) expression is unique to Gasoline. It is a member of the general family of pressure gradient expressions described in Monaghan (1992) (equation 3.5, setting σ = 1). Choosing this form was inspired by ideas presented in Ritchie & Thomas (2001) aimed at minimizing errors in strong density gradients. It has the desirable property that it closely mimics the form of equation (10), where the local density (1/ρ_i) can be taken out front and the two pressures have the same multipliers. It is found to minimize surface tension type effects highlighted by Agertz et al. (2007), as is demonstrated in the tests in Section 3. It also naturally complements the form of the energy equation.

The internal energy equation uses a form similar to that advocated by Evrard (1988) and Benz (1990),

\begin{eqnarray} \frac{{\rm d}\,u_i}{{\rm d}t}& = & \sum _{j} m_j \left( \frac{P_i}{\rho _i \rho _j} + \frac{1}{2} \Pi _{ij} \right) \boldsymbol {v}_{ij} \cdot \bar{\nabla _i W_{ij}}, \end{eqnarray}

(13)

where u_i is the internal energy of particle i and |$\boldsymbol {v}_{ij}=\boldsymbol {v}_i-\boldsymbol {v}_j$| is the velocity vector difference between particles i and j. Note the re-use of the Geometric Density Average. In concert with equation (12) this form conserves energy exactly in each pairwise exchange. For a purely adiabatic gas this form is equivalent to P/ρ multiplied by an SPH estimate for the divergence and thus closely follows the physical equation (11).

2.3.1 Entropy conservation

SPH models the fluid system via pairwise exchanges of momentum and energy between particles. These exchanges must be symmetric to achieve momentum, angular momentum and energy conservation. This conservation is only exact in the limit as the time-step approaches zero and is otherwise dependent on the quality of the integrator. Exact momentum conservation can be achieved with any integrator and fixed time-steps. If the SPH force expressions are derived from a Lagrangian (as in Springel, Di Matteo & Hernquist 2005; Hopkins 2013), symmetry is automatic and the form of the equations is fixed. This approach limits ones ability to adjust the form to minimize problems such as surface tension-like behaviour at sharp jumps in the pressure. On the other hand, the Lagrangian approach provides accurate entropy conservation Springel & Hernquist (2002). It is instructive to show how this comes about. In the absence of shocks, sources or sinks, the entropy function, A(S) = u/ρ^γ − 1 is constant on a mass element. For this to remain true, density and internal energy must evolve consistently so that

\begin{eqnarray} -\nabla \cdot \boldsymbol {v} = \frac{1}{(\gamma -1) u} \frac{{\rm d}u}{{\rm d}t} = \frac{1}{\rho } \frac{{\rm d} \rho }{{\rm d}t}. \end{eqnarray}

(14)

For an arbitrary SPH formulation, even if energy is manifestly conserved, there may be a systematic offset in the estimate for du/dt so that when density changes by orders of magnitude the thermal energy does not change appropriately and entropy is not conserved. One approach is to evolve entropy directly but this can create energy errors. The issue can be addressed directly by correcting the momentum and energy equations to exactly respect equation (14). In formulations derived from a Lagrangian, this correction, called f_i in Springel & Hernquist (2002), is automatically introduced by the Lagrange multipliers.

Our internal energy equation (equation 13) is only dependent on the local particle pressure making it closer to the physical energy equation (11) than most alternative forms. As a consequence entropy is closely conserved.

However, it is possible to explicitly estimate a correction for it (or any SPH formulation) by comparing the estimate for |$\nabla \cdot \boldsymbol {v}$| on the left-hand side of equation (14), arising from the internal energy (equation 13), to that on the right-hand side of equation (14) derived from equation (3),

\begin{eqnarray} \frac{1}{\rho _i} \frac{{\rm d}\,\rho _i}{{\rm d}t}& = & \sum _{j} \frac{m_j}{\rho _i} \frac{W^{\prime }}{h_i^5} \boldsymbol {v}_{ij} \cdot \boldsymbol {r}_{ij} \end{eqnarray}

(15)

\begin{eqnarray} & - & \sum _{j} \frac{m_j}{\rho _i} r_{ij}^2 \left[ \frac{\,W^{\prime }}{h_i^6} + \frac{3\, W}{h_i^4}\ \right] \frac{{\rm d}h_i}{{\rm d}t}. \end{eqnarray}

(16)

The term inside the square brackets in equation (16) comes from ∂ρ/∂h and was used to derive correction factors f in Springel & Hernquist (2002). Such correction terms are commonly referred to as ‘grad h’ terms. Under the idealized integral approximation of the sum, this term is identically zero. We find no systematic effect on entropy integration arising from it for the formulation presented here. In cases where severe errors occur (large amounts of expansion or contraction), the dominant effect is from the first term above. In equation (13), ρ_j replaces ρ_i in this term. This gives the form of f_i presented in equation (8).

Using the correction in equation (8) ensures that the ratio of the numerical estimates of the two key terms in equation (14) does not systematically vary from unity. Therefore entropy is closely conserved. This is demonstrated in the test section below, such as for the Evrard test where Gasoline exactly follows the correct entropy evolution prior to the shock and in the Sedov test. It can also be easily demonstrated by simulating an expanding sphere for many expansion factors and showing that the entropy is conserved. Traditional SPH forms without correction terms (e.g. the original Gasoline, Monaghan 1992) see systematic variations in the entropy.

2.4 Shocks and artificial viscosity

All hydrodynamics codes require extra numerical dissipation to broaden dissipative shocks whose natural width is smaller than their resolution element. Without this they are numerically unstable. This regime is common in astrophysics where molecular viscosity is very low. Codes based on piecewise reconstructions and Riemann solvers (e.g. Woodward & Colella 1984; Springel 2010b; Hopkins 2015) do this through slope limiters, making the solver first order in shocks and thus locally equivalent to the original Godunov (1959) scheme. Central schemes (e.g. Kurganov & Tadmor 2000) apply slope limiters but add dissipation in a manner similar to artificial viscosity. Direct artificial viscosity is the dominant approach in SPH. In all schemes, the goal is to minimize dissipation away from shocks. The original SPH form (Monaghan 1992) is as follows:

\begin{eqnarray} \Pi _{ij} = \left\lbrace {\begin{array}{l@{\quad}l}\frac{-\alpha \frac{1}{2}(c_i+c_j)\mu _{ij}+\beta \mu _{ij}^2}{\frac{1}{2} (\rho _{i}+\rho {j})} & {\rm for }\,\boldsymbol {v}_{ij}\cdot \boldsymbol {r}_{ij}< 0,\\ 0 & {\rm otherwise}, \end{array}}\right. \end{eqnarray}

(17)

where

\begin{eqnarray} \mu _{ij} = \frac{h(\boldsymbol {v}_{ij}\cdot \boldsymbol {r}_{ij})}{\boldsymbol {r}_{ij}^{\,2}+0.01 (h_i+h_j)^2}, \end{eqnarray}

(18)

where c_j is the sound speed of particle j. α = 1 and β = 2 are standard values for parameters representing shear and Von Neumann–Richtmyer (high Mach number) viscosities, respectively. We have implemented this and also the signal speed variant of the artificial viscosity (Monaghan 1997) in Gasoline2,

\begin{eqnarray} \mu _{ij} = \boldsymbol {v}_{ij}\cdot \hat{r}_{ij}. \end{eqnarray}

(19)

The viscosity is then effectively a function of the signal speed |$v_{\rm sig, ij} = \frac{c_i+c_j}{2} - \min (\boldsymbol {v}_{ij}\cdot \hat{r}_{ij}, 0),$| which better approximates a local wave-speed of the Euler equations.

In past work with Gasoline, we used the multiplicative Balsara (1995) limiter to reduce artificial viscosity,

\begin{eqnarray} \xi _{\rm Balsara} = \frac{\mid \nabla \cdot \boldsymbol {v}\mid }{\mid \nabla \cdot \boldsymbol {v}\mid +\mid \nabla \times \boldsymbol {v}\mid }. \end{eqnarray}

(20)

This limiter is effective in non-shocking, shearing environments but is otherwise very noisy and is thus half on (∼0.5) in most of the simulation volume.

Morris & Monaghan (1997) suggested a time-dependent limiter implemented via evolving the viscous parameter α, such that |${\rm d}\alpha /{\rm d}t = \max (0,- \nabla \cdot \boldsymbol {v}) - \alpha /\tau$|⁠. This increases in compressions and, otherwise, slowly decays with τ ∼ 0.1 h/c. The decay acts somewhat like a smoothing in time and helps compensate for noise. CD pointed out that α increases gradually, potentially too slowly for strong shocks. A second problem is that it cannot differentiate between compressions and shocks. CD suggested instantly setting α to a maximal value based on the time derivative of the divergence, |$-{\rm d}(\nabla \cdot \boldsymbol {v})/{\rm d}t$|⁠. This approach allows for rapid activation of artificial viscosity. It is noisy away from shocks, particularly when unavoidable particle reordering is occurring such as during density changes. This tends to make the limiter also partly on in much of the volume (e.g. α ∼ 0.15 between the rarefaction and the contact discontinuity in shock tubes in CD Fig. 11). Thus, the fact that CD chose not to apply a minimum α may be a moot point. As noted by Rosswog (2015) artificial viscosity in disordered regions could be of net benefit and provide a better overall solution if it helps limit velocity noise. For this reason, a noise-sensitive estimator can still perform well on many tests and Rosswog advocated for an explicit noise-based limiter.

A key problem with any viscosity approach based on divergence (including the CD approach) is that it creates viscosity in uniform compression (as seen in the left-hand panel of their fig. 15 where α ∼ 0.5 everywhere during the collapse test of Section 3.7). The usual indicator for shocks is |$\nabla \cdot \boldsymbol {v}$|⁠. For a spherical flow, |$\nabla \cdot \vec{v} = \frac{1}{r^2}\frac{\mathrm{\partial} r^2 v_r}{\mathrm{\partial} r} = \frac{\mathrm{\partial} v_r}{\mathrm{\partial} r} - 2 \frac{v_r}{r}$|⁠. In spherical collapse, divergence is commonly dominated by the second term and is not a reliable shock indicator. One could look for the most negative eigenvalue but in the collapse test, negative, non-radial eigenvalues dominate away from the shock.

A reliable indicator is the velocity gradient in the direction of the pressure gradient. This can be used in any geometry. We estimate the pressure gradient direction and velocity gradient using the local, one-sided estimator of Section 2.2 at the same time density is calculated,

\begin{eqnarray} \nabla P = (\gamma -1)\sum _j m_j u_j \nabla _i {\rm W}(r_{ij},h_i), \end{eqnarray}

(21)

\begin{eqnarray} \hat{n} = \frac{\nabla P}{|\nabla P|}, \end{eqnarray}

(22)

\begin{eqnarray} \frac{{\rm d}\,v}{{\rm d}n} = \sum _{\alpha ,\,\beta } n_{\alpha } V_{\alpha \beta } n_{\beta }. \end{eqnarray}

(23)

The result, |$\frac{{\rm d}\,v}{{\rm d}n}$|⁠, is an accurate local scalar shock indicator in the spherical collapse case but is still triggered in uniform compression. For this reason, we subtract off one-third the divergence (if it is negative) to ensure that the |$\hat{n}$| direction is playing a dominant role in the local compression. The resulting gradient based shock detector,

\begin{eqnarray} D = \frac{3}{2} \left(\frac{{\rm d}\,v}{{\rm d}n} + \max \left(-\frac{1}{3} \nabla \cdot \boldsymbol {v}, 0\right) \right). \end{eqnarray}

(24)

D takes on the extremal value of |$\frac{{\rm d}v}{{\rm d}n} = \nabla \cdot \boldsymbol {v}$| in the case of a shock but is only negative when the compression normal to the shock is more than one-third the overall compression. Other similar forms could be used. This indicator is also effective in the context of a Balsara-type limiter, replacing the divergence.

The general idea of limiters is to have a comparator that measures non-shocking flow gradients to combine with the shock indicator. Dehnen (private communication) noted that limiters such as that of Balsara (1995) that use the curl as a comparator lead to different shock properties in rotating systems. This motivated CD to use the magnitude of the trace-free shear tensor, |${\boldsymbol S}_{\alpha \beta }= 1/2 ({\boldsymbol V}_{\alpha \beta }+{\boldsymbol V}_{\beta \alpha }) - \delta _{\alpha \beta }\,\frac{1}{3} \nabla \cdot \boldsymbol {v}$| as the comparator. This is zero for pure rotation but still detects shear (e.g. differential rotation). Removing the trace makes it zero in isotropic compression or expansion that limits its ability to prevent viscosity there. In addition, in a strong shock the tensor is dominated by a single eigenvalue. If we align the shock normal with the x-axis then only the |${\boldsymbol S}_{xx}$| term is large and negative. Subtracting the trace makes the other two diagonal elements non-zero. Then, the norm |$|{\boldsymbol S}| = \sum _{\alpha ,\,\beta } {\boldsymbol S}_{\alpha \,\beta }^2$| is non-zero, which is undesirable. If we keep the trace and use |${\boldsymbol T}_{\alpha \beta }= 1/2 ({\boldsymbol V}_{\alpha \beta }+{\boldsymbol V}_{\beta \alpha })$|⁠, then an excellent overall indicator is |$\frac{D}{|{\boldsymbol T}|}$|

Our overall scheme to limit artificial viscosity is modelled on that of CD and is as follows:

\begin{eqnarray} \alpha _{{\rm loc},\,i} = \alpha _{\rm max} \frac{A_i}{A_i + v_{{\rm sig},i}^2}, \end{eqnarray}

(25)

where

\begin{eqnarray} A_i = 2\, h_i^2\,\xi _i \max \left(-\frac{{\rm d}D}{{\rm d}t}, 0\right). \end{eqnarray}

(26)

The only change here is in the definition of A_i, which is constructed with the 1D velocity derivative of equation (24). This combined with the different h definition in CD led us to boost the A_i term by 2. For the tests shown in this work, we employ α_max = 2. Experiments with strong shocks indicate higher peak α values may be needed. Our parameter values were chosen based on the test results shown in Section 3. We also note that we use the gradient estimator of equation (6) rather than the full treatment in the appendix of CD.

Whenever α is less than α_loc, it is set directly to that value; otherwise,

\begin{eqnarray} \frac{{\rm d}\,\alpha _i}{{\rm d}t} = (\alpha _{{\rm loc},\,i}-\alpha _i)/\tau _i, {\rm \ where\, } \tau _i = \frac{h_i}{0.2 c_i}. \end{eqnarray}

(27)

This decay is a little faster than that of CD. We have experimented with minimum α values but find that noise often contributes sufficient base α value where required. A more in depth examination is presented in Section 3.4 in the discussion of the Gresho–Chan test.

For the limiter we use,

\begin{eqnarray} \xi _i = \left(\frac{1-R_i}{2}\right)^4 \end{eqnarray}

(28)

where

\begin{eqnarray} R_i = \frac{\sum _j m_j D_j/|{\boldsymbol T}|_j W_{R,\,ij}}{\sum _j m_j W_{R,\,ij}}. \end{eqnarray}

(29)

R_i is constrained to be in the range [−1, 1]. Formally |$-1 \le \frac{{\rm d}v}{{\rm d}n}/|{\boldsymbol T}| \le 1$| by construction. However, D is a modification to |$\frac{{\rm d}v}{{\rm d}n}$| and noise occasionally pushes the value outside this range. The weighting W_R, ij could use the regular kernel except that this only strongly weights a small number of central particles and makes the estimate noisy. We use W_R, ij = (1 − (r_ij/(2 h_i))⁴. ξ_i is zero in expansions, ∼1/16 in intermediate and noisy regions (including uniform compression) and is maximized at 1 in shocks.

2.5 Diffusion

Diffusion for any scalar quantity in SPH, A_i, can be estimated using (Monaghan 1992),

\begin{eqnarray} \frac{{\rm d}A_i}{{\rm d}t}|_{\rm Diff} = -\sum _j m_j \frac{(d_i+d_j)(A_i-A_j)(\boldsymbol {r}_{ij}\cdot \nabla _i {\bar{W}_{ij}})}{\frac{1}{2}(\rho _i+\rho _j)\,{\boldsymbol {r}}_{ij}^2}, \end{eqnarray}

(30)

where d_i is the diffusion coefficient of particle i. If diffusion rates vary rapidly in space, a harmonic mean may be useful to replace (d_i + d_j)/2 above, but this is not our default.

As noted in Wadsley et al. (2008), traditional SPH lacks a mechanism to model the diffusion that is intrinsic to unsteady flows. Mesh codes diffuse numerically at rates proportional to the absolute fluid velocity that make them non-Galilean invariant (Wadsley et al. 2008; Hopkins 2013).

Turbulent diffusion dominates in high Reynolds number flows where shear can drive local fluid instabilities. Sub-grid turbulence may be modelled in SPH using the trace-free local shear tensor, |${\boldsymbol S}$| (Wadsley et al. 2008; Shen et al. 2010). We estimate the turbulent diffusion coefficient using |$d_i=C |{\boldsymbol S}| h_i^2$| with a coefficient C ∼ 0.03–0.1. We used C = 0.03 here. This estimate avoids diffusion in pure rotation and non-shearing flows. We use the shear tensor estimate described previously.

Every fluid scalar should diffuse, including thermal energy, u_i, metals and so forth, and this is done in Gasoline2. Diffusion should be considered an essential component of any SPH code unless the target is solely laminar flows. Incorrect results for entropy transport (e.g. galaxy cluster entropy profiles) and hydrodynamic instabilities are a consequence of its absence. We demonstrate this explicitly through the destruction of a gas blob via fluid instabilities in Section 3.3.

There are other astrophysically relevant sources of diffusion such as Spitzer conduction in hot plasmas (used in Keller et al. 2014). Monaghan (1987) has shown that artificial diffusion proportional to the viscosity can be necessary in extremely strong shocks. We not use it for the tests shown here. We apply all sources of diffusion within the same framework by adding to the diffusion rate coefficient, d_i.

2.6 Time integration

Gasoline has used Kick-Drift-Kick second-order leapfrog integration with hierarchical powers of two time stepping since the original implementation (Stadel 2001; Wadsley et al. 2004; Springel 2005). This approach is symplectic only if time-steps never change. When particles change time-steps (by factors of two each time), secular drift in otherwise conserved quantities can be introduced though this is not systematic and the code conserves energy well overall as shown in the original paper and confirmed in the tests. Recent work (Springel 2010b; Hopkins 2015) has suggested summing the total fluxes of every pairwise momentum and energy exchange to allow for exact conservation of these quantities. This is appealing, in principle, but not symplectic. Issues of this nature can be minimized by ensuring neighbouring particles have similar time-steps. Thus in Gasoline2, time-step criteria are applied in a pairwise fashion. For example, the Courant condition is applied as

\begin{eqnarray} \Delta t|_{\rm Courant, i} = 0.4 \min _j \frac{\bar{h}_{ij}}{1.25\bar{c}_{ij}+0.75\left[\alpha _{ij}\bar{c}_{ij} + \beta \mu _{ij}\right]}, \end{eqnarray}

(31)

which is symmetric between the particles. The term in square brackets arises where artificial viscosity is on. This form is standard (e.g. Monaghan 1992). We also employ the acceleration criterion |$\Delta t|_{\rm Acc} = 0.2 \sqrt{\frac{h}{a}}$|⁠.

Saitoh & Makino (2009) showed that large differences in time-steps between nearby particle can cause disastrous energy non-conservation. Following their suggestion, we never let a particle time-step exceed four times the current time-step estimate of any neighbour. In addition, if high speed events insert new neighbours with significantly shorter time-steps, we wake up the long time-step neighbours and cut short their current time-step to adjust it downwards as much as necessary. This makes the scheme temporarily first order but closely limits the overall error. It is critical for strong shocks such as the Sedov–Taylor blast, which has negligible initial temperatures (see Section 3.6).

We integrate heating, cooling and chemical networks using sub-cycling (arbitrarily many sub-steps, smaller than the hydro time-step). These equations are commonly stiff with very stringent stability requirements. Most work with Gasoline to date used solvers based on a modified semi-implicit Bulirsch–Stoer algorithm (Press et al. 1992) or the CHEMEQ2 scheme (Mott & Oran 2001). Algorithms based on code from Press et al. (1992) cannot legally be included in a public release. Gasoline2 also offers the option to cool using the public Grackle package (Bryan et al. 2014; Kim et al. 2014). These integrators use mid-point estimates of the hydrodynamic contributions so that the scheme is second order overall. We apply a time-step limiter to the hydrodynamics based on overall changes to the internal energy, |$\Delta t|_{\rm u} = 0.25\ u/(\frac{{\rm d}u}{{\rm d}t})$|⁠. The form of the energy equation (13) ensures that energies can never go negative due to errors in PdV estimates.

Grid codes, in particular, commonly use operator splitting for non-hydrodynamical terms (due to the modular nature of Riemann-solver approaches) and are thus first order in the energy integration when cooling is included because the cooling time-scale is typically much shorter. Such choices can have important implications for heating and cooling in astrophysical contexts and are relatively unexplored.

3 TEST PROBLEMS

In this section, we present a range of test problems on which we have run Gasoline2 and Changa. All results shown were made with Gasoline2. These tests have been selected because they are standard, they test the key problems exhibited by traditional SPH and show how the new methods resolve them. For our testing, we have used public initial conditions where available. However, we also insist on using three dimensions, glass-type initial conditions and one set of standard parameters with all modern SPH features included. This includes optimizations such as multiple time-steps which can negatively impact conservation of energy.

Unlike grid codes, tests in one or two dimensions do not accurately reflect how SPH will perform in 3D simulations. 1D tests avoid interpenetration issues and re-gridding noise. 2D simulations can be set up in a stable close-packed hexagonal initial grid, which minimizes noise. A relaxed glass configuration is the natural state for structure that is spontaneously formed in simulations, such as through collapse or instabilities. Glasses always have some noise and cannot be completely static.

For any numerical method, it is possible to massage the results by changing code parameters for different tests. A well-known example of this is slope limiters, which can be pushed into an antidiffusive extreme to avoid spreading at contact discontinuities (Woodward & Colella 1984). However, this choice causes numerical instabilities in production simulations where a more conservative choice is preferred. Thus, in the following tests, we sometimes show how different components of the method affect the results, but we always conclude with a demonstration of how the full Gasoline2 method performs with standard parameters. We have also avoided excessive resolution so the tests measure practical performance. We show the results warts and all, so to speak.

3.1 Square test

This test is used to demonstrate the presence of surface tension-like effects (Hopkins 2013; Saitoh & Makino 2013; Hu et al. 2014; Hopkins 2015). We perform this test in three dimensions starting with glass initial conditions similar to the set-up of Hu et al. (2014): 256³ equal-mass particles in a unit cell with a square region of four times higher density in −0.25 < x < 0.25 and −0.25 < y < 0.25. 2D versions and regular lattices favour un-representative results for this test particularly. Fig. 1 shows a slice through the result at time t = 3, with a thickness equal to the particle spacing for Gasoline2 in its standard configuration and a version re-run with traditional force terms. The traditional SPH version clearly suffers from a surface tension-like problem. The key method improvement that gives this result in the use of the GDFs (Section 2.3). These results are quite comparable to those of Hu et al. (2014) and demonstrate that pressure–entropy formulations (Hopkins 2013) are not required to alleviate surface tension. Note that the initial condition is not in perfect force balance at the corners, which results in a small correction (rounding) there. Saitoh & Makino (2013) demonstrated that resolving this issue completely requires variable particle masses and uniform volume elements per particle, essentially giving up density-based adaptivity.

Figure 1.

A comparison between Gasoline2 (which uses Geometric Density Forces), left-hand side, and a run with tradition SPH forces, right-hand side, at time t = 3 on the square test. The red particles are those initially within the square region.

3.2 Kelvin–Helmholtz instability

The Kelvin–Helmholtz instability is a direct test of how surface tension and insufficient mixing suppress mixing via fluid instabilities. We use the initial condition based on that used in Read et al. (2010),² We note that those initial conditions feature particles aligned on a regular lattice, rather than a glass, so we generate a set of initial conditions with identical resolution and properties using three glass slabs. Aside from the glassy initial particle positions, our initial conditions are identical to those used by Read et al. (2010). The domain of this IC is a rectangular slab (L_x, L_y, L_z) = (1, 1, 1/32). This slab is composed of three regions, a central slab where |y| < 0.25, and two slabs above and below this. The volume is initialized with uniform pressure, where ρ_central/ρ_outer = T_outer/T_central = 2. These flows are shearing relative to each other, with each moving in opposite directions with v_x = 0.11 c_s. A sinusoidal perturbation in v_y is imposed on each of the two boundaries, defined by equation (32) where δv_y = 4v_x and λ = 0.5.

\begin{eqnarray} v_y = \delta v_y [\sin {(2\pi (x+\lambda /2)/\lambda ))}\exp {(-(10(y-64))^2)} \nonumber\\ &&- \,\sin {(2\pi x/\lambda )}\exp {(-(10(y+64))^2)}] \end{eqnarray}

(32)

We generate our ICs using a set of three glasses, rather than a grid or lattice, to avoid artificially suppressing noise introduced by the natural SPH re-gridding.

As shown in Fig. 2, using GDFs is sufficient to eliminate surface tension and allow the instability to grow. The addition of diffusion produces growth rates in better agreement with high-resolution simulations. The full code version shown in the right-hand panels includes the variable viscosity limiter. As can be seen in the figure, these are not essential for this test. We note that alternate, stronger forms of diffusion such as the artificial conduction of Price (2008) also improve results on this test even with more traditional pressure gradients. However, we would argue that modified gradients with relatively low turbulent diffusion is more generally useful.

Figure 2.

Kelvin–Helmholtz instability test. The most dramatic improvement comes from the use of GDFs. The left-hand panel shows results with just this change from traditional SPH. The central panel includes diffusion, and the right-hand panel is the full Gasoline 2 code result.

3.3 Blob test

One of the primary difficulties traditional SPH faces is mixing of multiphase fluids. The ‘Blob’ test (Agertz et al. 2007) is designed to study this, by embedding a stationary dense, spherical cloud in a uniform supersonic flow. The cloud itself is in pressure equilibrium with the flow, and should be broken up by the development Kelvin–Helmholtz and Rayleigh–Taylor instabilities as it is accelerated up to the flow velocity.

The initial conditions for the Blob test come from the Wengen 3 test suite. The domain is a periodic rectangular prism, with dimensions (L_x, L_y, L_z) = (10, 10, 40) in units of the cloud radius, which is centred at (x, y, z) = (0, 0, −15). This cloud is placed in pressure equilibrium, with ρ_cloud/ρ_wind = 10. The wind velocity is thus |$v_{{\rm wind}} = \sqrt{10}R_{{\rm cloud}}/\tau _{{\rm KH}}$|⁠.

The evolution of a density slice with the full Gasoline2 code is shown in Fig. 3. The blob develops fluid instabilities, breaks up and diffuses into the flow as expected. Agertz et al. (2007) quantified this result by examining the fraction of the original blob mass still above 64 per cent of its original density as shown in Fig. 4. The upper panel of this figure shows the progression with a linear vertical axis (as in the original paper) for different versions of the code. The lower panel shows the progress on log axes. This plot shows how traditional SPH suffers with surface tension, which inhibits the initial break up. Adding grid-scale diffusion modestly improves the result.

Figure 3.

Agertz (2007) Blob Test: density evolution for the full Gasoline2 code (log grey-scale, factor of 200 in density).

Figure 4.

Agertz (2007) Blob Test: mass remaining versus time with comparison of methods. The top panel uses a linear vertical axis (as in presented in the original paper) and the lower is the same data on a log axis. The log y-axis emphasizes how traditional SPH keeps dense material for long periods while mesh codes and Gasoline2 break the material up at an accelerating rate.

On linear axes it looks like the GDF is sufficient to achieve good results. It follows the initial break up well. This period until time 3 or so is characterized by an accelerating break up as the pieces of the cloud break-up faster because their individual KH times are shorter. Thus on a log plot, the characteristic time to halve mass decreases and the slope becomes more negative. However, once the final pieces are close to the resolution scale (near time 3), a new issue arises that must be resolved through diffusion. The full Gasoline2 and older Gasoline results with the new force and diffusion are slightly different and inhabit a region of solution space similar to the range spanned by the enzo results. We have verified that other strong mixing models, such as the artificial conduction of entropy used in phantom (Price et al. 2017) produce similar results.

The enzo results include three curves because fixed-grid codes are not Galilean invariant. The v_blob = 0 curve is the standard run. In v_blob = 1, the blob moves to the left-hand side and the background flow is static and for v_blob = 0.5 each had half the motion. SPH gives identical results in these three cases. However, for grid codes there is faster breakup when the flow moves rapidly relative to the grid. This result emphasizes that there is no unique answer to this problem. It is sensitive to initial conditions, method choices and effective resolution. Though the initial condition was symmetric, all the codes develop asymmetries as the instabilities go non-linear.

In Fig. 5, we show entropy slices in a zoomed-in regions around the blob for several different methods on the test at different times. The top case is the standard enzo result. The blob test initially has uniform high entropy in the flow and low entropy in the blob. As the flow impacts the blob the shock creates even higher entropy in the flow, shown in white in the figure. The blob is somewhat like a mushroom cloud on its side, where deceleration plays the role of gravity and the high entropy gas initially wraps around the exterior. The shock rapidly fades and from that point new entropy production is minimal. In traditional SPH, the lack of diffusion means that SPH particles retain their entropy values for the rest of the run. This is particularly apparent in the second and fourth rows where no diffusion was used. In particular, even with surface-tension removed (GDF), there are dense, low entropy knots that move through the high–low entropy material. Resolution-scale instabilities rapidly mix entropy on a local crossing-time. The precise numerical behaviour depends on the absolute velocity in grid codes (three enzo cases) and the mixing model (in SPH). The primary difference between the second last and final rows is the move to more neighbours (64 to 200) and the Wendland C4 kernel. The result is a reduction in noise and smoother features.

Figure 5.

Blob Test results showing entropy for different methods (T^3/2/ρ on log colour scale, factor of 1000 range, bright is high). The top row shows standard enzo (AMR-PPM) results. The second row is Gasoline (2004) and the third adds just diffusion. The fourth row shows just the additions of the GDF and the fourth shows that with diffusion. The bottom row shows results with the full Gasoline2 method. These results correspond to the curves shown in Fig. 4.

3.4 Gresho–Chan vortex

The Gresho–Chan vortex (Gresho & Chan 1990) is a cylindrical vortex that is ideal for testing how well a code can simulate an inviscid fluid. The ideal solution is steady, with centrifugal acceleration balancing pressure gradient. Viscosity will cause angular momentum transport, disrupting this equilibrium by letting the inner part of the vortex torque up the outer parts. The vortex is initialized with ρ = 1, −1 < x, y < 1, and a piecewise pressure function,

\begin{equation} P(r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}5+12.5\,r^2 & \quad (0 \le r < 0.2)\\ 9+12.5\,r^2-20\,r+4\ln {5r} & \quad (0.2 \le r < 0.4)\\ 3+4\ln {2} & \quad (r \ge 0.4)\\ \end{array}\right. \end{equation}

(33)

and velocity function, v_r = 0, and

\begin{equation} v_\theta (r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}5r & \quad (0 \le r < 0.2)\\ 2-5r & \quad (0.2 \le r < 0.4)\\ 0 & \quad (r \ge 0.4).\\ \end{array}\right. \end{equation}

(34)

We evolve the vortex to t = 3, or ∼2.4 rotations of the peak. We note in passing that not all code papers do this (e.g. Hu et al. 2014). Given that the vortex decays in all cases, it is mostly important for making fair comparisons. A key choice was to evolve this problem using a glass in three-dimensions, where most others (Springel 2010b; Kawata et al. 2013; Hopkins 2015) evolve the problem in two. We refer the reader to Dehnen & Aly (2012) for a detailed discussion of this problem and how numerical factors affect it.

Fig. 6 shows the results at four different resolutions at time t = 3. A viscosity limiter (Morris & Monaghan 1997, CD) is critical for this test. The improvements present in the Gasoline2 limiter (Section 2.4) do not change the outcome from a more straightforward CD limiter on this problem. The convergence of the L1 norm error is fairly typical and somewhat sub-linear (⁠|$N_{1D}^{-0.8}$|⁠) as shown in the left-hand panel in Fig. 7.

Figure 6.

Velocity versus radius in the Gresho–Chan Vortex with the full Gasoline2 method at t = 3 (2.4 rotations of the peak). Resolution is indicated in the top left-hand side of each panel. The blue bars are inserted at one particle separation as an indication of resolution and the thickness is an indicator of the rms deviation. The grey points are individual particles. The black curve is the exact solution.

Figure 7.

Gresho–Chan Vortex L1 errors. Left-hand panel: L1 error versus effective resolution for N_1D = 32, 64, 128 and 256. The convergence rate is roughly as |$N_{{\rm 1D}}^{-0.8}$|⁠. The Noise result is nearly on-top of the base result. Right-hand panel: L1 errors versus amount of artificial viscosity introduced to combat noise for 0.1 < r < 0.3 (near the peak of velocity). The solid lines show the total L1 error and the dashed curves show the contribution to L1 error from the binned (mean) solution only.

Rosswog (2015) argued that triggering viscosity on noise could improve the solution further. We experimented with noise triggers similar to those of Rosswog (2015) as additions to the shock trigger of CD. For all the variants, we found that it was hard to achieve a substantial improvement over our standard scheme. We include results of additional tests with a very simple noise trigger, |$\alpha _{{\rm NOISE}}=\sigma _v^2/(\sigma _v^2+c^2)$| where σ_v is the rms velocity noise at the particle to illustrate the general behaviour (light line in the left-hand panel in Fig. 7). Attempting strong noise suppression made the solution worse (diamond line, 10 α_NOISE). We attribute this partly to the fact that the shock trigger also triggers on noise in the absence of strong expansion or shear to suppress it so some noise viscosity is already present. Examining our results, we were also concerned that it was possible to reduce the L1 error even though the mean solution (i.e. in radial bins) was worse. This behaviour is demonstrated in the right-hand panel in Fig. 7, which measures the error near the peak of the vortex. We chose the peak because there are many particles outside r = 0.3 that otherwise dominate the measurement. The L1 norm is a combination of the spread and the mean deviation from the exact solution. We find that extra viscosity always makes the mean solution worse (squares in the right-hand panel in Fig. 7) even if the overall L1 error is modestly improved. The velocity distributions with the noise trigger are visually indistinguishable from the results shown in Fig. 6. For now the standard Gasoline2 has no noise trigger but a careful re-examination might yield benefits from employing one.

3.5 Sod shock tube

No code paper would be complete without the classic Sod shock tube (Sod 1978). This simple 1D Riemann problem begins with two domains at rest, but out of pressure equilibrium. The initial conditions are ρ_left = 1, P_left = 1 for x < 0 and ρ_right = 0.25, P_right = 0.1795 for x > 0. The self-similar solution consists of a rarefaction fan travelling to the left-hand side, a contact discontinuity and a shock moving to the right-hand side and is shown in the figures as a black line.

Fig. 8 shows Gasoline2 without variable viscosity but with the Balsara limiter. This is comparable to standard Gasoline as of 2011 with the addition of the Wendland C4 kernel and 200 neighbours, α = 2 and β = 4. The use of 200 neighbours and the large shock parameters increases the shock width over versions with older kernels, 50–64 neighbours, α = 1 and β = 2. Choosing α = 1 substantially decreases the shock width and allows some modest post-shock ringing, as seen in Wadsley et al. (2004).

Figure 8.

Sod shock tube for Gasoline2 without variable viscosity. The symbols show averages values separated by the local 1D particle spacing to indicate the resolution and the thin lines indicate the exact solution. The top panel shows density and the lower is velocity.

Fig. 9 shows the full Gasoline2 result, including the viscosity limiter on same initial conditions as Fig. 8. The solution to the left-hand side of the initial contact has very little viscosity and experiences mild ringing at the velocity peak of the rarefaction wave. All variants of the time-dependent viscosity suffer in the rarefaction wave where there is no viscosity until after it peaks. Applying a minimum alpha was not found to be helpful in reducing this unless it was large enough to negatively affect other tests (e.g. Gresho–Chan Vortex). The method supplies an α ∼ 0.1–0.2 in the post-rarefaction region, as was also seen by CD on this test. We note that the velocity in this region for PSPH also shows some mild ringing (Hu et al. 2014; Hopkins 2015).

Figure 9.

Same Sod shock tube as Fig. 8 for Gasoline2 with the viscosity limiter. The additional lower panel shows the local viscosity α.

3.6 Sedov–Taylor Blast

The Sedov–Taylor blast wave (Taylor 1950; Sedov 1959) is a spherical shockwave generated by the injection of energy in a central region. For this test, we have a 128³ box, with 64 particles injected with 6.78 × 10⁵³ erg. The background fluid has temperature, T = 0, making this an infinite Mach number shock, with a density enhancement of 4 for a γ = 5/3 equation of state. The initial conditions have a domain of 6 × 6 × 6 kpc, and a density of 0.5 cm^{− 3}. This test will break codes without careful time-step adjustment, as demonstrated in Saitoh & Makino (2009). The total energy error between the start and the final time in the Gasoline2 run was 0.2 per cent.

Fig. 10 shows results for Gasoline2 at three times. For the final time, we also show the velocity and entropy function profiles. The initial condition had a finite initial hot mass with a maximum entropy indicated by the diamond in the right-hand panel so the entropy follows the exact solution until it gets to this value. In the left-hand panel, we plot the density on log axes at three equally spaced times. The solution deviates near the point where the entropy limit kicks in. The central panel shows velocity that is commonly left out in plots of these results. As shown in Hu et al. (2014), we see considerable ringing in the velocity. We have verified that the viscosity rises to α = 4, its maximum value, well before the shock. The perfectly cold (T = 0) pre-shock gas essentially forces this. This is an infinite Mach number shock and thus a strenuous test of the code. We have experimented with larger α_max values but they both broaden the shock and negatively affect other aspects of the solution for this and other tests. Thus, we show the current results as an acceptable compromise.

Figure 10.

Sedov–Taylor explosion results for zero external temperatures. This is both an extremely strong shock and a test for accurate time integration. The panels show density at three times (left), velocity at the final time (centre) and entropic function at the final time (right). Note that the exact solution (black lines) is for a point explosion.

3.7 Evrard collapse test

The adiabatic gas-only collapse test from Evrard (1988) includes a strong shock and cold pre-shock infall. It is the only test shown here that requires gravity. It can be used to assess the ability of a code to follow the collapse phase during the formation of an astrophysical object. It tests shock capturing and the coupling between gravity and gas dynamics. We use a glass initial condition with 28 000 particles, comparable to the resolution in the original papers. The error in total energy over the course of the full Gasoline2 run was 0.04 per cent

Fig. 11 shows the structure at time t = 0.8 (t = 0.88 in Hernquist & Katz 1989). The results shown as diamonds are binned at the particle spacing with actual particle values shown as lighter points. The solid line is a high resolution 1D PPM solution from Steinmetz & Mueller (1993). Prior SPH results suffer from pre-shock entropy production as the compression triggers the artificial viscosity. In SPH codes that we are aware of (including Gasoline in 2004), roughly half the entropy is produced during the pre-shock infall because viscosity is triggered on negative divergence.

Figure 11.

Evrard adiabatic collapse test at time t = 0.8 with standard Gasoline2 parameters. The four panels show density, velocity, entropy and properties relevant to shock capturing. Note that entropy does not increase until the particles are within range of the actual shock near r = 0.1.

Using a gradient based shock indicator based on |$\frac{{\rm d}v}{{\rm d}n}$| rather than divergence, we detect the shock at the right place. The extra shock width is solely due to the finite resolution. The difference is highlighted in the top right-hand panel of Fig. 11, where we show that our shock viscosity α turns on at the right place, where |$\frac{{\rm d}v}{{\rm d}r}$| switches sign from positive to negative. As discussed in Section 2.4, in a collapse like this the divergence is dominated by the geometric contraction of the solid angle rather than the radial velocity gradient. Using gradient-based shock detection works well in this and any geometry. It also filters out uniform contraction or expansion.

Fig. 12 shows Gasoline2 results with less conservative smoothing (Wendland C2 kernel and 64 neighbours). This improves the resolution (central density and narrower shock) at the cost of increased noise everywhere, particularly in velocity. This illustrates that in large-scale collapse problems with orders of magnitude increases in density, the benefits of smaller local errors associated with large neighbour counts are less obvious.

Figure 12.

Evrard adiabatic collapse test at time t = 0.8 except with the Wendland C2 kernel and 64 neighbours.

4 CONCLUSIONS

This paper provides a complete description of the SPH methods present in Gasoline2. We argue that these improvements or similar ones are required for any modern SPH code to overcome the limitations of traditional SPH. In the following summary we have highlighted methods unique to Gasoline2 with italics:

• GDF and internal energy expressions which both minimize surface tension effects in multiphase flow and naturally provide excellent entropy conservation.

• Turbulent Diffusion, first introduced in Wadsley et al. (2008), which models physically necessary sub-grid turbulent mixing, alleviating problems with hydrodynamic instabilities and scalar transport (e.g. entropy).

• A Gradient-Based shock detector based on the 1D velocity gradient normal to shocks, rather than divergence. This prevents spurious viscosity in convergent flows such as spherical collapse problems.

• Time-dependent artificial viscosity, modified from the prescription of Morris & Monaghan (1997) and CD, to reduce velocity noise and better treat inviscid flow.

• Wendland kernels, shown by Dehnen & Aly (2012) to allow larger neighbour numbers without pairing and improve convergence.

• Time-step adjustment, shown by Saitoh & Makino (2009) to be necessary for handling sudden jumps in time-step such as due to rapid heating in strong shocks.

We have presented test problems that show the necessity of each of these components to give satisfactory results. The tests presented here cover many situations relevant to astrophysical simulations. Thus Gasoline2 and Changa should be regarded as cutting-edge, modern SPH codes. It is worth emphasizing that we use a Density–Energy SPH method. Pressure–entropy formulations are not required prevent surface tension problems. In fact, the distinction between these method has been blurred via the use of many different density weightings (e.g. Rosswog 2015).

Lagrangian techniques such as SPH offer considerable advantages for astrophysical simulations (e.g. natural adaptivity, Galilean invariance, accurate orbits and efficient time-stepping). We note that SPH is still somewhat more diffusive and slower to converge than other methods for a given resolution. New methods have arisen that offer faster convergence through the use of more complex gradient estimates and fluxes from Riemann solvers, such as arepo (Springel 2010b) and gizmo (Hopkins 2015). arepo suffered from some SPH-like numerical convergence issues, which have now been resolved (Mocz et al. 2015; Pakmor et al. 2016). New SPH methods with accurate gradients based on integral form (ISPH) also show considerable promise without the cost of using Riemann solvers (Cabezon, Garcia-Senz & Figueira 2016). Valdarnini (2016) showed promising ISPH results for previously problematic regimes such as subsonic turbulence. Thus, continuous code improvement and testing new methods seem to be the lot of the computational astrophysicist for the foreseeable future.

Acknowledgements

The authors would like to thank Hugh Couchman, Walter Dehnen, Jonathan Panuelos and Zachariah Levine for helpful interactions during the preparation of this manuscript. We would also like the thank the referee, Daniel Price, for helpful suggestions. The authors would also like to thank SHARCNET/Compute Canada for HPC resources and support and the Natural Sciences and Engineering Research Council of Canada for financial support. Tom Quinn acknowledges support from NSF award AST-1311956.

1

Gasoline2 is a public code, and is available for download at http://gasoline-code.com

2

That IC, and those for the blob test, are available for download at http://astrosim.net/code/doku.php?id=home:codetest:hydrotest