Improving smoothed particle hydrodynamics with an integral approach to calculating gradients
^{1}
Dept. de Física i Enginyeria Nuclear, Universitat Politècnica de Catalunya,
Compte d’Urgell 187,
08036
Barcelona,
Spain
^{2}
Institut d’Estudis Espacials de Catalunya, Gran Capità 24,
08034
Barcelona,
Spain
^{3}
Departement Physik. Universität Basel, Klingelbergstrasse 82,
4056
Basel,
Switzerland
email: ruben.cabezon@unibas.ch
Received:
23
August
2011
Accepted:
14
November
2011
Context. The smoothed particle hydrodynamics (SPH) technique is a wellknown numerical method that has been applied to simulate the evolution of a wide variety of systems. Modern astrophysical applications of the method rely on the Lagrangian formulation of fluid Euler equations, which is fully conservative. A different scheme, based on a matrix approach to the SPH equations is currently being used in computational fluid dynamics. These matrix formulations achieve better interpolations of the physical magnitudes but they are, in general, not fully conservative. The matrix approach to the Euler equations has never been used in astrophysics.
Aims. We develop and test a fully conservative SPH scheme based on a tensor formulation that can be applied to simulate astrophysical systems.
Methods. In the proposed scheme, derivatives are calculated from an integral expression that leads to a tensor (instead of a vectorial) estimation of gradients and reduces to the standard formulation in the continuum limit. The new formulation improves the interpolation of physical magnitudes, leading to a set of conservative equations that resembles those of standard SPH. The resulting scheme is verified using a variety of wellknown tests, all of them simulated in two dimensions. We also discuss an application of the proposed tensor method to astrophysics by simulating the stability of a Sunlike polytrope calculated in three dimensions.
Results. The proposed scheme is able to improve the results of standard SPH in the twodimensional tests, especially in the simulation of subsonic hydrodynamic instabilities. Our results for the stability of the Sunlike polytrope suggest that the new method can be used in astrophysics to carry out threedimensional calculations with a computational cost that is only slightly higher (i.e. ≤50% for a serial code) than that of a standard SPH formulation.
Conclusions. A formalism based on a matrix approach to Euler SPH equations was developed and checked. The new scheme is more accurate because of the renormalization imposed on the interpolations, which is fully conservative and probably less prone to undergo the tensile instability. The analysis of several test cases suggest that the method may improve the simulation of both subsonic and supersonic systems. An application of the tensor method to astrophysics is, for the first time, successfully carried out. These encouraging results indicates that more work should be invested in the applications of matrix SPH formulations to astrophysics.
Key words: methods: numerical / hydrodynamics / instabilities
© ESO, 2012
1. Introduction
The hydrodynamic method known as smoothed particle hydrodynamics (SPH) is a gridless Lagrangian approach to continuum mechanics devised by Gingold & Monaghan (1977) and Lucy (1977). As demonstrated in an endless amount of applications to computational fluid dynamics, ranging from largescale astrophysical simulations to nuclear physics on the Fermi level, SPH is a robust and confident way to simulate the dynamical evolution of deformed systems. The success of SPH relies strongly on the way gradients are estimated. For example, the value of density in a given spatial coordinate is obtained from particles located in the neighborhood using an interpolating function called kernel. As the kernel is an analytical differentiable function, it is easy to obtain gradients just evaluating the gradient of this weighting kernel function. Then it is quite straightforward to write the Euler equations of fluid mechanics in terms of the kernel and its derivatives (Monaghan 1992, 2005).
Despite the success, SPH has also a number of shortcomings and weak points that make the technique less useful than it should be to handle a number of wellidentified problems. First of all, the numerical noise is usually larger than in other techniques, which poses a problem for applications involving very subsonic movements, as in the case of hydrodynamical instabilities. Another difficulty has to do with the use of the artificial viscosity (AV) formalism to handle shock waves and smooth the inherent presence of noise. Several solutions and recipes have been invoked to handle these problems, none of them totally satisfactory, among them: variable AV coefficients (Morris 1997), artificial heat fluxes to smooth pressure discontinuities (Price 2008), nonstandard approximations to the momentum equation (Abel 2011), and a refined treatment of kernel interpolations as in the technique known as moving least squares interpolation (i.e. MLS; Dilts 1999).
In this paper, we present an approach to the momentum and energy equations that combines a novel way to estimating gradients and the variational principle, and allows an improved modeling of physical phenomena, especially in systems subjected to small perturbations. The proposed formalism, based on an integral approach to the derivatives (IAD), includes the standard formulation as a particular case. It leads to a matrix formulation similar to that of MLS methods, but the form of the SPH equations in IAD remains closer to that of standard SPH. The text is organized as follows. In Sect. 2 the mathematical formalism for calculating gradients is described and discussed. In Sect. 3 we develop a formulation of the momentum equation compatible with the variational EulerLagrange principle in the SPH framework, and propose a consistent and fully conservative energy equation. Section 4 is devoted to validating the hydrocode through the simulation of four known tests where standard SPH has, to a smaller or larger degree, difficulties. In Sect. 5 we elaborate an extension of the method that by combining the exact estimation of the derivative of linear functions and good momentum and energy conservation, can be applied to simulate linearized systems. The ability of the proposed method to simulate the structure of real threedimensional (3D) astronomical bodies is discussed in Sect. 6. Finally, the main conclusions of our work, as well as some comments about the shortcomings of the developed scheme and future lines of improvement are outlined in the last section, which is devoted to conclusions.
2. Integral approach to first derivatives
There is an ample literature describing the treatment of the first and second derivatives in SPH (Monaghan 2005; Rosswog 2009). In the following we take a different route to estimating first derivatives which leads to a more general expression than the most common extant derivative procedures. Our method resembles to those known as MLS methods where interpolations relative to gradient estimation were constrained to enhance linearity by imposing a minimization procedure on the errors.
Inspired by the treatment of second derivatives (Brookshaw 1985) we define an integral expression from which we will estimate gradients. As in the case of second derivatives, it is expected that an integral approach will lead to a less noisy estimation of first derivatives. Therefore, we define (1)where f(r) is any differentiable function and W(r′ − r,h) is a spherically symmetric interpolating kernel, usually a sharplike Gaussian of width h (where h is usually referred as the smoothing length). The size of h is often interpreted as the local resolution of the simulation. Expanding the factor [f(r′) − f(r)] to first order (2)and writing (3)Eq. (1) can be expressed as a matrix equation: (4)where (5)and (6)Once the linear system in Eq. (4) is solved we get the value of the gradient of the function f. In the case of spherically symmetric kernels, τ_{11} = τ_{22} = τ_{33} and τ_{ij} = τ_{ji} = 0, ij. Taking for instance the Gaussian kernel (7)where n is the dimensionality of space and , it leads to . Therefore, for the Gaussian kernel the gradient of f reduces to (8)which is the same expression as Eq. (2.13) of Monaghan (2005) even though it has been obtained through a different procedure. Since the standard kernels are spherically symmetric Eq. (8) can be simplified to (9)which in SPH parlance becomes (10)
which because of its simplicity and good performance has become the most popular procedure for calculating first derivatives in SPH. Therefore, Eq. (10) is a particular case of IAD, Eq. (4).
Nevertheless, many of the above symmetry properties of tensor are lost when integrals are converted into finite summations. For example, the elements on the diagonal of matrix in Eq. (4) do not necessarily have the same value and elements outside the diagonal can differ from zero. However, the formulation of the first derivative in terms of the matrix Eq. (4) has a very interesting feature: the derivative of a linear function is always exact by construction. The demonstration of such property is straightforward in 1D SPH: (11)taking f_{a} = px_{a} + q,f_{b} = px_{b} + q the expression above gives .
In two dimensions, the explicit solution for Eq. (4) is written as (12)where (13)Unfortunately, direct application of the tensor scheme above leads to SPH equations of movement that are not fully conservative. Nevertheless, a slight modification of the scheme makes it possible to build exact conservative equations by taking advantage of the symmetry properties of the kernel to discretize the vector I(r) as (14)Thereafter, the formulation that results after substituting Eq. (14) into Eq. (4) is labeled IAD_{0}. This restricted interpretation of Eq. (4) implies that exact evaluation of the derivative of linear functions cannot be achieved if we are to also preserve the exact conservation of momentum and energy. Nevertheless, we present below strong indications that using IAD_{0} leads to a better evaluation of the derivative of linear functions than the standard procedure given by Eq. (10), especially when a small or moderate number of neighbors is used to compute summations. We can illustrate this point in the following numerical experiments.
As a first test, we simulate a static system where we compare the gradient of a linear distribution of density in one dimension obtained with IAD, IAD_{0} and standard schemes as a function of the number of neighbors n_{b}. A linear density profile ρ(x) = 1 + x was obtained aligning N = 100 equidistant particles of adequate mass along the xaxis. The density was calculated using (15)Note that even though IAD should provide the precise derivative, in practice it never does owing to the small errors in the calculation of density. It is also worth noting that the standard calculation of the derivative has two additional potential sources of error. For a linear density profile ρ_{b} = ρ_{a} + p(x_{b} − x_{a}), and using Eq. (10), we have (16)where p is the “exact” derivative and ϵ_{1} ≃ 0, ϵ_{2} ≃ 1. On the other hand, there is only one source of error, ϵ_{3} ≃ 0, if IAD_{0} is used, given by (17)In Fig. 1 (upper and bottom left) we represent the value of the relative error in the derivative , with respect to the analytical value p as a function of the smoothing length h, normalized to the interparticle distance Δ, for the Gaussian and cubic spline kernels respectively. Independently of the kernel, the error when full IAD linear interpolation was used is always much smaller than for the other methods. When IAD_{0} or the standard derivative were used, the error increased appreciably. The error is large when the smoothing length is shorter than the interparticle distance but decreases rapidly when h increases, as expected. Nevertheless, the error in the standard calculation always remained larger, especially for a small or moderate number of neighbors (n_{b} ≤ 5) regardless of the kernel we used, but especially for the cubic spline. Since the particle sample is highly ordered, the error curve of the Gaussian kernel is much smoother than that of the cubic spline owing to the infinite range of the kernel. We verified that when the Gaussian is truncated to 2h, the error profile becomes qualitatively similar to that of the cubic spline. However, the error profile for IAD_{0} follows almost exactly the error curve for the gradient of density calculated as a simple quotient for adjacent particles. In this sense, IAD_{0} seems to be more coherent with the computed density distribution of the sample than IAD, a trend that also holds in 2D, as shown below.
Fig. 1 Relative error ϵ in the first derivative of density ρ(x) = (1 + x) g cm^{3} as a function of the smoothing length (in interparticle units) calculated using both tensor and standard SPH schemes. Crosses are for direct derivative estimation. Upper left is for the error when the Gaussian kernel is used. The upper right picture shows the contribution of several error sources to the total error ϵ. Bottom pictures are the same but for the cubicspline kernel. 

Open with DEXTER 
The contributions to the total error of ϵ_{1},ϵ_{2} and ϵ_{3} in Eqs. (16) and (17) are shown in the rightmost part of Fig. 1. The contribution of ϵ_{1} and ϵ_{3} are similar but the error ϵ_{2}, which affects only the standard scheme and exhibits a large oscillation in the case of compact supported interpolators, represents the main channel to the total error in the derivative. In general, tensor methods become more reliable for higher dimensionality. Thus, we conducted a similar numerical experiment in more than one dimension where the matrix approach can potentially be more beneficial.
As a second test, a bidimensional system was set using a sample of N = 62500 equally spaced particles on a twodimensional lattice. The mass of the particles along the xaxis was conveniently modified to reproduce a linear density profile, ρ(x,y) = 1 + x. The first derivative of the density, for different number of neighbors was obtained using Eqs. (4) and (10). In the upper row of Fig. 2 we show the relative error in the derivative, , with respect to the analytical value p as a function of the smoothing length. Calculations were carried out using the Gaussian and the cubic spline kernels. The trend here is similar to that of the 1D calculation but magnitudes fluctuate less, especially in the case of the cubic spline. When the number of neighbors is large the error is small in all three methods, but tensor schemes are more accurate than the standard procedure. Nevertheless, for shortening smoothing lengths, h ≤ 1.5Δ, the error curve steepens and for h < Δ standard SPH begins to give unacceptable results while the tensor estimations of the derivative remain applicable. This result also agrees with the 1D calculation shown in Fig. 1, which shows that working with n_{b} ≤ 3,4 in 1D may give an appreciable error in the derivative of a linear function if Eq. (10) is used to compute the derivatives. As in the 1D case, the gradient calculated as a direct fraction (crosses in Fig. 2) remained closer to the results obtained using IAD_{0}.
Fig. 2 Relative error ϵ in the first derivative of linear functions calculated in 2D as a function of the smoothing length. Upper rows depict the value of ϵ for ρ(x,y) = (1 + x) g cm^{3} obtained using the Gaussian (left) and the cubic spline (right) kernels. Bottom rows are the same but for the function g cm^{3}. Error profiles with crosses are for the direct derivative estimation. 

Open with DEXTER 
Finally, we considered a spherically symmetric linear profile ρ(x,y) = 1 + r, where . In the bottom row of Fig. 2, the value of the relative error ϵ is shown as a function of h(Δ). The profile of the error follows a slightly different trend than in the previous case. As before, the error is large when the smoothing length is shorter than the interparticle distance but in the case of the Gaussian kernel there is a minimum in the error profile at h ≃ 1.1Δ. From that point on the error smoothly increases. This behavior is not so clear in the case of the cubic spline where the error curve roughly stabilizes above h = 2Δ. There is, therefore, an optimal value of h in this case that minimizes the errors in the first derivative. Close to the center linearity is lost because of the symmetry imposed on the distribution. When the value of h is large, interpolations are affected by that geometrical constraint and, at some point, the advantages of working using a large smoothing length are lost. Nevertheless, the errors calculated using IAD_{0} are always smaller than that of the standard scheme. Again the gradient of density estimated as a simply quotient is more closely reproduced by IAD_{0} than by the other schemes.
3. The momentum and energy SPH equations using IAD_{0}
Nowadays the most accurate formulation of the momentum equation in SPH comes from the variational principle. It has been shown (Monaghan 2005, and references therein) that the solution of the EulerLagrange equations leads naturally to a very symmetric scheme including the effects of spatial gradients on the smoothing length. The resulting equation conserves momentum by construction and, equally important, ensures perfect energy conservation for nondissipative systems. It is shown below how the tensor approach built using the IAD_{0} scheme can also be compatible with the variational principle.
The Lagrangian of the system is evaluated as (18)where v_{b}, u_{b} and s_{b} are the velocity, specific internal energy and specific entropy of particle b. Using this Lagrangian and assuming isentropic evolution, the EulerLagrange equations subjected to the constraint ρh^{3} = const., lead to the equation of movement (Springel & Hernquist 2002)^{1}(19)where is a term accounting for the gradient of h. The icomponent in Eq. (19) can be written (20)An estimation of the density gradients using the tensor Eq. (4) is (21)where . Elements of matrix are those defined in Eq. (5) after changing integrals to summations. For particle a, they are defined to be (22)and the jcomponent of vector I in the IAD_{0} approach is (23)Thus, for particle a the icomponent of density gradient is (24)where c_{ij,a} are the elements of matrix associated with particle a, and d is the dimension of the space. Explicit expressions for these elements in cartesian 2D were given by Eq. (13). Equation (24) can be used to directly compute the second term on the right of Eq. (20).
We can evaluate, in the same way, for particle b(25)To calculate we use Eq. (15): (26)The last term in Eq. (26) also appears during the evaluation of ∇_{b}ρ_{b}, namely (27)Comparing Eqs. (27), (26) and (25), we finally get (28)Substituting Eqs. (24) and (28) into Eq. (20), the icomponent of the momentum equation for particle a is given by (29)where It is then straightforward to derive the appropriate energy equation (32)Since the magnitudes , defined by Eqs. (30) and (31), are antisymmetric, the conservation properties of the IAD_{0} formulation are identical to those of standard SPH.
As in standard SPH, an artificial viscosity term has to be added to stabilize the numerical scheme and handle shocks. The AV increases the effective pressure in those zones of the fluid that undergo compression. Including the viscous acceleration, the momentum equation is given by (33)where Π_{ab} accounts for the viscous pressure (34)where the symbols have their usual meaning and r_{ab} means r_{a} − r_{b}. The coefficient μ_{ab} is (35)To preserve momentum conservation, we evaluate the arithmetic mean of (36)The inclusion of artificial viscosity in the energy equation leads to (37)As a representative example of the ability of the new scheme to handle subsonic physics, we consider the stability of an inhomogeneous system in mechanical equilibrium. A sample of N = 62500 particles were evenly distributed within a square lattice with periodic boundary conditions ensuring that the initial density was set to ρ_{0} = 1 g cm^{3}. The density was then perturbed at random allowing maximum percentage variations of 5% across the lattice. The internal energy of each particle was adjusted to ensure isobaricity with P_{0} = 1 dyn cm^{2} leading to a sound speed value c_{s} ≃ 1.3 cm s^{1}. The evolution of the system was followed during almost a sound crossing time, t_{s} ≃ 0.5 s using the three schemes, IAD, IAD_{0} and standard SPH, with n_{b} = 30 and 100. Calculations using full IAD were carried out using the Eqs. (52) and (53) described in Sect. 5. Hereafter, the word standard, STD, refers to the modern, fully conservative, Lagrangian formulation of SPH that explicitly includes the gradient of the smoothing length parameter in the scheme (see for instance Rosswog (2009) and references therein). The results of the simulations are summarized in Fig. 3, where we present the evolution of the standard deviation in the pressure P(t) with respect to P_{0} for the three aforementioned schemes. For this particular test it is clear that both tensor schemes provide a more persistent mechanical equilibrium than standard SPH during a sound crossing time, especially when full IAD is used. As expected, the standard deviation in the pressure decreases as the number of neighbors increases.
Fig. 3 Evolution of magnitude of an inhomogeneous 2D system, initially in hydrostatic equilibrium, calculated using the different SPH schemes mentioned in the text. 

Open with DEXTER 
3.1. Computational issues
The increase of the computational requirements of the integral method is more evident now. With respect to memory requirements, the nine coefficients (in 3D) of matrix and those belonging to matrix can share the same memory space. Thus, nine coefficients per particle have to be stored to compute in Eq. (33). While the increase in memory requirements are not a significant problem, the burden in the computational time is larger. First of all, one has to compute six of the nine coefficients of symmetric matrix , Eq. (22). Although these calculations involve simple operations, they have to be performed in a separate treewalk because previous knowledge of density is necessary. Afterwards, matrix has to be inverted in order to calculate the coefficients c_{ij,a} for particle a. Fortunately, these matrix inversions do not take much time because they can be completed out of the tree structure. Finally, one has to compute the nine quantities of coefficients associated with each pair of neighbors particles a,b; instead of the three that are necessary in the standard formulation, and use them to compute both the momentum and the energy equations.
It is obvious that the proposed method demands a larger computational effort than the standard SPH. Nevertheless, the extra burden could be small if the physical problem under simulation requires the calculation of longrange forces (i.e. gravity) or involves complex physics (i.e chemical or nuclear reactions, transport phenomena). In addition note that, according to Figs. 1 and 2, for standard SPH to achieve a similar accuracy as IAD_{0} in calculating the first derivative, the number of neighbors has to be larger. To maintain the spatial resolution, the total number of particles also has to be larger by the same factor. Therefore, changing both the number of neighbors and total number of particles could have a larger impact on the computational performance than using IAD_{0}.
4. Hydrodynamic tests
The aim of this section is to check the theory described in Sects. 2 and 3 using four tests, all of them widely used in CFD. The first two simulations concern the growth of the KelvinHelmholtz (KH) and RayleighTaylor (RT) hydrodynamic instabilities, whose evolution always remains subsonic and, as we demonstrate, the fully conservative IAD_{0} scheme leads to improved results with respect to the standard method for the same identical initial conditions. The last two tests are related to strong supersonic flows where shock waves take over, as in the wall heating shock and Sedov tests, where the use of the tensor route to estimate derivatives again provides more accurate results, especially for the wall heating shock test.
The simulations described in this section were carried out in a cartesian two dimensional scenario because of the higher achieved resolution, easiest initial setting and analysis of the results. Periodic boundary conditions were imposed in all tests except in the RT case, where particles from then top and the bottom of the box were fixed. A perfect gas equation of state (EOS) P = (γ − 1)ρu, where u is the specific internal energy and γ = 5/3, is always used. The coefficients α, β of artificial viscosity, in Eq. (34), were set to 1 and 2 respectively, except for the wall heating test where we took α = 1.5, β = 3. The cubic spline was used to perform the SPH interpolations but several calculations carried out with other interpolators (i.e. harmonic kernels with index n = 3; Cabezón et al. (2008)) did not lead to significant differences. We used a selfadaptive smoothing length parameter which keeps constant the number of neighbors, n_{b}, of a given particle. Although the simulations described below were obtained for n_{b} = 100, the same calculations carried out using n_{b} = 36 neighbors did not provide significant differences. Corrections due to the gradient of the smoothing length were explicitly included in the IAD_{0} formulation as it comes from the LagrangeEuler variational principle.
In Table 1 we show a summary of energy and momentum conservation properties for the different calculated models. The tensor calculations are always more able to ensure momentum and energy conservation than the standard ones for the same elapsed time. The lack of energy conservation during the pointlike Sedov explosion in all methods is due to the hard initial conditions. In contrast, the very good conservation of momentum in the supersonic numerical experiments is partially due to the spherical symmetry imposed to the initial configurations.
Conservation properties of the computed models: KelvinHelmholtz (KH), RayleighTaylor (RT), the wall heating shock (WH) and Sedov (SED).
4.1. KelvinHelmholtz instability
The KelvinHelmholtz instability is a wellknown test to check the ability of any hydrodynamic code to handle subsonic perturbations. This instability appears when there is a sufficient velocity shear in the interface layer between two fluids with different densities. Small perturbations of the velocity field in the orthogonal direction to the interface grow, leading to a mixing of both fluids. This is usually simulated in a box with periodic boundary conditions where two fluid regions are defined with densities ρ_{1} and ρ_{2} respectively. Both layers have opposite parallel velocities leading to a shear discontinuity in the contact interface. To develop the instability, a small perturbation is seeded in the interface as a sinusoidal mode of length scale λ.
In our case, we have simulated a central band of high density fluid (ρ_{1}) moving in a lowdensity medium (ρ_{2}) in a squared lattice of 1 cm side in the XY plane using N = 62500 particles. The mass of the particles was arranged to obtain the correct density profile following a ramp function Robertson et al. (2010). In this way, we smoothed the interface density jump to make it comparable to the SPH resolution using (38)where A is a normalization constant and Δy = 0.05 cm. Therefore, the density profile was given by (39)where ρ_{1} = 2 g cm^{3} and ρ_{2} = 1 g cm^{3}.
The seed of the perturbation is obtained using a sinusoidal function for the v_{y} component of the velocity field. For the initial velocity, we then have (40)where we assume n = 2 and Δv_{y} = 0.1 cm s^{1}. Also note that the v_{x} component has been smoothed using the same ramp function used for the density, and that v_{1} = 0.5 cm s^{1} and v_{2} = −0.5 cm s^{1}, which corresponds to the high and low density bands respectively.
Figure 4 shows four snapshots of the growth of the KelvinHelmholtz instability at different times for the calculation using IAD_{0} (first row) and the standard SPH implementation (second row). As can be seen, the standard formulation does a poor job of resolving the structure of the instability. Although in the beginning the main shape is in gross agreement with the tensor simulations, the final image is blurry, the interface is not welldefined and the shape is incorrect. In the case of the tensor calculation, the instability grows cleanly and at good rate, and the definition of the extremes of the billows, where the finest structure appears, is clearly enhanced. To achieve similar results to those obtained with the IAD_{0} technique, different methods have been proposed to maintain the standard description, mainly based on including an artificial thermal conductivity (Price 2008). This method provides reliable results, but includes a new set of parameters and estimates, such as maximum signal velocity between two particles to obtain a diffusion parameter. Furthermore, the inclusion of an artificial thermal conductivity leads to an extra dissipation of gradients away from discontinuities, hence it needs to design some means of controlling the amount of dissipation.
Fig. 4 Evolution of the KelvinHelmholtz instability. The snapshots show the color map of density at times t = 0.1,1.0,2.0 and 3.0 s for methods IAD_{0} (first row) and standard SPH (second row), with Δv_{y} = 0.1 cm s^{1}. Times t = 1.0,3.0,4.0 and 5.0 s for methods IAD_{0} (third row) and standard SPH (fourth row), with Δv_{y} = 0.01 cm s^{1}. 

Open with DEXTER 
In Fig. 5, we present the evolution of total linear momentum during the development of the instability. Although both schemes give a very satisfactory momentum conservation, we can see that conservation using the IAD_{0} method is at least one order of magnitude better than that of standard scheme.
To test the tensor approach in a more hostile scenario, we diminished the value of the amplitude of the initial velocity perturbation by an order of magnitude (i.e. Δv_{y} = 0.01 cm s^{1}). In the third and fourth rows of Fig. 4 we show the results of the simulations using the configurations IAD_{0} and STD. It is clear that in the standard formulation the instability was unable to grow, while it does using IAD_{0}. Changing the geometry of the initial particle distribution from a square to a uniform hexagonal lattice or reducing the size of the smoothing length taking n_{b} = 36 neighbors did not appreciably change the results. We can then state, that the matrix approach has inherent virtues that are absent in the standard formulation, mainly related to the generalization of the derivation technique to a tensor expression, which, in some sense, diminishes the errors derived from the discretization.
It is also worth to mention that during the simulations carried out using the standard calculation a generalized clumping of particles in the low density regions often appeared. This is a wellknown problem, called pairing or tensileinstability, where particles get “stuck” if owing to an unstable stressstrain relation that occurs typically when the second derivative of the kernel becomes negative and high strain is developed between particles. This problem eventually leads to difficulties in resolving the distances between neighbors, incorrect interpolations, and finally large nonphysical behavior of the particles that halts the calculation. A typical way to cope with this problem is to adaptively change the kernel slope so that it becomes more centrally peaked (Cabezón et al. 2008). However, the IAD_{0} method was unaffected by this problem, hence we found that the tensor method seems to be less susceptible to pairing instability. The ability of matrix methods to cope with pairing instability has been reported by other authors (Oger et al. 2007), and they merit further investigation in the context of IAD_{0}.
Fig. 5 Evolution of the absolute deviation from initial linear momentum during the development of the KelvinHelmholtz instability. Figure on the left is for the x component of total momentum calculated using IAD_{0} and standard (STD) schemes, whereas figure on the right is the same for the y component. 

Open with DEXTER 
4.2. RayleighTaylor instability
The simulation of the growth of the RayleighTaylor instability in a stratified fluid in the presence of gravity is a classical test of subsonic fluid dynamics. In its simplest form, a box containing two fluids with different densities separated by a sharp transition region is placed inside a gravitational field (Youngs 1984). If the denser fluid is located on top of the lighter, the system becomes physically unstable and any small perturbation of the interface leads to fluid overturn. The denser fluid sinks and the lighter one rises, producing characteristic structures known as spikes and bubbles respectively. After a linear phase of growth, which can be studied analytically, the system enters a complex nonlinear regime. In real fluids, however, the viscosity of the fluid prevents the growth of the instability at short wavelengths, thus viscosity tends to dampen or even to completely suppress the instability. One serious problem of numerical schemes incorporating the artificial viscosity formulation is that they introduce too much viscosity into the system, compromising the growth of instabilities. As we later show, the integral approach to the derivative, Eqs. (33) and (37), provides a substantial improvement in the simulation of the RT phenomenon when the initial perturbation amplitude is small.
A sample of N = 62500 particles was distributed in a squared lattice of 1 cm side. The mass of the particles was conveniently arranged to reproduce the density profile given by (41)where y is the vertical coordinate with the origin at the interface and Δy is the width of the transition zone between both fluids. The density in the two zones were set to ρ_{1} = 1 g cm^{3} and ρ_{2} = 2 g cm^{3} respectively, and the size of the transition zone was taken to be Δy = 0.05 cm.
The integration of the hydrostatic equilibrium equation gives the pressure distribution along the fluid (42)where g is the gravity, here of value g = 0.5 cm s^{2}. With this setting, the speed of sound at the interphase is c_{s} = 0.7 cm s^{1}. A perturbation in the boundary layer with λ = 0.25 cm was seeded by giving an initial small vertical velocity v_{y} to the particles according to the prescription (Abel 2011) (43)and v_{y} was set to zero for y positions above 0.7 cm and below 0.3 cm. The value of velocity perturbation was set to Δv_{y} = 0.01 cm s^{1}. In Fig. 6 we show several snapshots of the development of the Rayleigh Taylor instability calculated using IAD_{0}. As we can see, the RT instability is able to grow after a few tenths of a second. When the same calculation was attempted using the standard SPH scheme to compute the acceleration, the results were of much lower quality, the instability was totally damped and there was no development of the RT fingers. As in the case of the KelvinHelmholtz instability, changing the initial particle setting from the square to an hexagonal regular lattice or reducing the size of the smoothing length did not qualitatively alter the above conclusions. A similar conclusion was reported by Abel (2011) who used a variant of standard SPH, albeit not fully conservative, to simulate the grow of the RT instability using a similar initial setting of the experiment. This negative result does not mean that standard SPH is unable to simulate this phenomena because it could handle it after a careful initial setting that minimizes the numerical noise. Nevertheless, the results of our simulations indicate that the tensor method is less dependent on the particular geometry of the initial lattice.
Fig. 6 Development of the RayleighTaylor instability for Atwood number 1/3 calculated using the IAD_{0} scheme. The snapshots show the color map of density at times t = 0.4,4.2,5.1 and 5.7 s. 

Open with DEXTER 
In Fig. 7, we show the evolution of the center of mass of the spikes y_{s} from the initial contact surface. According to the theory, the evolution is exponential during the initial stage ξ(t) = ξ_{0}exp[Γ(t − t_{0})], roughly until the initial perturbation ξ_{0} at t_{0} (here we took t_{0} = 1.1 s) has grown to a size ξ_{l} ≃ 1/k where k is the wave number. In our tests, ξ_{l} ≃ 0.04 cm. The theoretical growth rate Γ_{t} during the exponential phase is given by , where A_{t} = (ρ_{2} − ρ_{1})/(ρ_{2} + ρ_{1}) is the Atwood number. The exponential phase is followed by a brief second stage in which the penetration of the dense fluid evolves as y_{s} ∝ gt^{2} until drag takes over and both falling spikes and rising bubbles reach a limiting speed. According to Fig. 7, all these features are present in the simulations. Nevertheless, there was only a qualitative agreement with the analytical growth rate during the initial exponential phase because the value of Γ inferred from the simulations was almost one half of the theoretical value. This is unsurprising because a reliable quantitative estimation during the initial stage needs both a higher resolution as well as a better initial setting than that we are currently using. In particular, the transition zone between light and dense fluids has to be sharper to correctly represent the Atwood number in this stage. An analytical expression for the terminal velocity of a rising bubble in a tube was obtained by Layzer (1955), where l is the radius of the tube. Taking l = λ = 0.25 cm, we get v_{b} = 0.104 cm s^{1} in good agreement to the numerical value deduced from Fig. 7.
Fig. 7 Evolution of the center of mass of the spikes during the growth of the RayleighTaylor instability depicted in Fig. 6. The velocity of the center of mass of the bubbles, v_{b}, and spikes, v_{s}, is also shown. For t > 4 s, a limiting speed is reached. 

Open with DEXTER 
4.3. The wall heating shock test
We now consider simulations of supersonic events. In the test known as the wall shock problem, a spherical or cylindric supersonic stream of gas is launched towards its geometrical center forming a highly compressed region. For these geometries, simulations can be compared to the predictions of an analytical approach to the evolution of thermodynamic variables as a function of the initial conditions (Noh 1987). Although the gross features of the event are correctly captured by SPH simulations, it is wellknown that schemes based on artificial viscosity have difficulties in providing a detailed description of the wall heating shock test. The reason is that AV spreads the shock over several computational cells, which induces a nonphysical increase in the internal energy ahead of the shock. For converging flows, a large artificial spike in internal energy is observed at the convergence center. As a consequence, a profound dip in the density profile appears to keep the pressure smooth. However, the inclusion of a good amount of artificial viscosity is mandatory to providing a successful description of the shock.
Our purpose behind this test was to discern whether the integral route to calculating derivatives can provide more reliable results than the standard procedure. Therefore, we use IAD_{0} and Eqs. (15), (33) and (37) to carry out the simulations.
A sample of N = 57600 particles of the same mass were evenly spread across a lattice to ensure that the density was homogeneous, ρ = 1 g cm^{3}. The initial pressure and internal energy of particles were negligible. The system was imploded by imposing a spherically symmetric velocity field v_{r} = −1 cm s^{1}. In Fig. 8 we show the density and radial velocity profiles at t = 0.3 s, when the shock is welldeveloped. Both profiles follow a simple pattern that remains close to the analytical estimations. In the central region, r < 0.1 cm, a plateau of compressed material forms, with a characteristic density of ρ_{s} ≃ 14.5 g cm^{3} a little lower than the analytical value of ρ = 16 g cm^{3} for γ = 5/3 (Noh 1987). In this inner region matter remains stagnated with a velocity close to zero. Outside the plateau the density abruptly declines through the shock front trying to regain its initial value, whereas the velocity increases to reach v_{r} = −1 cm s^{1} at the incoming, but still unshocked, matter. As we can see in Fig. 8, the simulation using the tensor approach is of higher quality than that of standard derivatives. In particular, the numerical oscillations in the postshock region in both, density and velocity, were considerably smaller in the tensor formulation. Nevertheless, the dip in the density profile and the maximum value achieved by density were similar in both calculations.
Fig. 8 Density and velocity profiles during the wall heating shock test. Figure on the left is for IAD_{0} calculation, whereas figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 16 g cm^{3}. 

Open with DEXTER 
4.4. Sedov test
In the Sedov test, the evolution of a shock wave front born as a consequence of a pointlike explosion is studied as it propagates in a homogeneous medium. The problem of an intense explosion in a gas is a standard test for hydrocodes that is of relevance to astrophysics, where it is not rare to find strong shocks in many scenarios involving fluid motions at high velocity. The theoretical solution was found by Sedov by applying selfsimilar methods and dimensional analysis for different geometries (Sedov 1959). In its simplest formulation, the Sedov problem has an initially cold gas at rest. At t = 0 s there is a point explosion at the origin that in Sedov (1959) was treated as an instantaneous release of energy at the origin, and assumed that the background material through which the expanding gas sweeps behaves like a perfect gas. For these initial conditions there are precise, albeit algebraically complicated, analytical expressions for the fluid variables. In the case of SPH, the artificial viscosity smears the shock over 2–3 times the smoothinglength. As a consequence, the density jump across the shock front is always smaller than the factor of four predicted by the theory for γ = 5/3. In more than one dimension, resolution issues are crucial not only to resolve the peak of the blast wave but also to reproduce the correct postshock variables downstream and the structure of the rarefied tail close to the origin. We want to investigate the extent to which the use of the integral approach to the derivative can improve the results of the SPH simulations.
Fig. 9 Density and velocity profiles during the Sedov test. Figure on the left is for IAD_{0} calculation, whereas the figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 4 g cm^{3}. 

Open with DEXTER 
The initial setting was similar to that of the wall heating test but this time the velocity was set to zero everywhere at t = 0 s, to ensure that the initial configuration is in equilibrium. To create the explosion, a large amount of internal energy was instantaneously released in a small region around the center of the box. To smooth the initial discontinuity, we took an initial pressure step that decays as a Gaussian function (44)where P_{1} and P_{2} are the pressures in zones left and right of the step and σ sets the width of the pressure decay. In the present simulations, P_{1} = 10^{4} dyn cm^{2} and P_{2} = 1 dyn cm^{2} and σ^{2} = 36 cm^{2}, which smooths the internal energy step over about 5–6 times the smoothing length. In Fig. 9 we depict the density and velocity profiles once the selfsimilar state has been achieved. As we can see the profiles calculated using the tensor and standard approaches to calculating gradients do not differ so much. The peak in the density profile is about 80% of what is expected from a strong shock moving through a γ = 5/3 gas but the profile in the postshock zone is smoother in the calculation using IAD_{0}. The velocity profile is also similar but the tensor calculation showed fewer oscillations and a more ordered behavior in the postshock tail. These results agree with our main conclusions given in the previous test dealing with the wall heating shock, reinforcing the idea that for identical initial settings IAD_{0} provides more reliable results than standard derivatives also for supersonic events.
5. Handling linear phenomena: momentum and energy equations using IAD
Ensuring both exact linear interpolation and perfect momentum and energy conservation using either full IAD or other similar schemes, seems to be difficult. In this respect, the best physical systems to try to conciliate both things are pure linear systems for which linear acoustic wave propagation is a representative example. A suitable alternative SPH formulation can be found starting from the differential acceleration equation (45)which, using Eq. (10), can be easily accommodated within the SPH formalism (46)This equation can be deduced independently using the variational principle and, provided the kernel is symmetrized , it conserves linear and angular momentum exactly. If instead of Eq. (10) we use Eq. (8) to compute Eq. (45), the resulting equation is (47)This particular form of the acceleration equation has several interesting features: 1) In the continuum limit it reduces to the standard expression, Eq. (46), because the last term in the equation vanishes. 2) As shown in Sect. 3.2, an isobaric system with smooth density gradients remains in equilibrium with negligible acceleration for a long time (see Fig. 3).
However, it is also evident that the inclusion of the last term in Eq. (47) breaks its symmetry. Still, in the linear limit, momentum is conserved to a high extent as can be easily demonstrated in 1D. We first expand the magnitudes P / ρ^{2} in Eq. (47) assuming a smooth density gradient (48)where Δρ = ρ_{a} − ρ_{b} and Δρ / ρ ≪ 1. Inserting Eq. (48) into Eq. (47) and neglecting higher order terms, the acceleration is given by (49)If we assume a linear pressure profile within the kernel domain of the particle (50)the acceleration becomes (51)which is the onedimensional Newton equation, thus leading to an exact value for the acceleration. In consequence, global momentum is wellpreserved provided that the system remains in the linear regime and the density gradients are smooth. According to this discussion, we complete the acceleration equation obtained in Sect. 3 using the IAD_{0} scheme with a corrective term (52)where 0 ≤ Ψ ≤ 2 is a parameter controlling the strength of the applied correction. For Ψ = 0, the already checked IAD_{0} scheme results, whereas for Ψ = 2 we have full linear IAD.
To take into account the corrective term in the energy equation we include the power released/absorbed by such term in Eq. (37) (53)such that total energy is exactly conserved. For very subsonic movements, the energetic contribution of the correction term is small. For Ψ = 0, Eqs. (52) and (53) reduce to the fully conservative IAD_{0} scheme. Alternatively, the case of Ψ = 2 could be used to handle linear systems with smooth density gradients as it shows enhanced interpolation abilities and good conservative properties.
Sound wave propagation is an example of a system that evolves in the linear limit of Euler equations. It is therefore a test ideally suited to highlighting the potential advantages of the IAD scheme. A homogeneous system with ρ_{0} = 1 g cm^{3} was simulated using N = 62500 particles uniformly distributed in a lattice of size 1 cm. The initial values of both pressure and density were set to one. A sample of n_{p} = 200 particles inside a circle of radius R_{p} = 0.032 cm located at the center of mass of the lattice was obliged to oscillate with small amplitude and period P. This compact set of particles emulated the effect of an external oscillating piston moving adiabatically onto the initially unperturbed gas. As the piston moves, it launches regular trains of circular waves that propagate at the sound speed cm s^{1} for γ = 5/3. We want to check the ability of the different aforementioned schemes to handle with this problem.
Particles belonging to the piston move homologously, following the harmonic oscillator law (54)where ω stands for the angular frequency, and the maximum value of the radial velocity of a particle located at distance r from the center is (55)The resolution places limits on both the period and the maximum velocity at the piston head . The value of the period was set to P = 0.05 s so that around eight complete waves were launched before the first one reached the limits of the system. The value of the maximum velocity at the piston head was set to cm s^{1}.
Fig. 10 Acoustic wave profiles at t = 0.35 s calculated using IAD, IAD_{0} and standard SPH (STD). The waves were generated by the periodic displacement of a circular piston of size 0.035 cm located at the center of the lattice. The continuum line is the analytical solution calculated taking a sound speed of c_{s} = 1.29 cm s^{1}. 

Open with DEXTER 
The results of the simulations for the three SPH approaches: IAD, IAD_{0} and STD are summarized in Fig. 10, which represents the radial velocity profile of the gas at time t = 0.35 s. It is quite evident that the spherical symmetry was better preserved during the IAD calculation, while there is a progressive degradation of the symmetry in the IAD_{0} and standard calculations respectively. The worst case corresponds to the standard scheme but even in this case spherical symmetry is reasonably preserved. A more quantitative analysis can be done making use of the analytical solution for diverging circular waves, usually known as waves in a membrane in the literature. The solution for free harmonic traveling waves in a circular membrane is written as (56)where u = kr, k is the wave number and J(u) is the Bessel function with radial and angular modes m = n = 0. The Bessel functions are constrained by J(0) = 1 at the center of the wave train. The profile for the radial velocity given by Eq. (56) with and the comparison with the numerical profiles obtained using the different schemes are also shown in Fig. 10. The analytical solution reminds one of a damped cosine wave propagating at the sound speed c_{s} = 1.29 cm^{1}. On the whole, the three schemes were able to correctly track the evolution of the waves. However, a close inspection reveals small differences among them. The best value of wave velocity corresponds to the IAD_{0} scheme ( cm s^{1}) followed by the standard ( cm s^{1}), whereas full IAD falls a bit short ( cm s^{1}). Even though the damping in the three schemes is always larger than that of the analytical solution, it seems that the worst case also corresponds to the fully IAD scheme built setting Ψ = 2 in Eqs. (52) and (53).
The simulation of linearized systems represents the most favored situation to apply exact linear interpolation schemes such as IAD because, as shown above, conservation of momentum and energy is probably as good as in the standard formulation. Nonetheless, the numerical experiments with the acoustic waves suggested limited advantages of using these schemes. In fact, some details of the wave evolution are better represented by IAD_{0} or even by the standard scheme. Only the spherical symmetry was more accurately preserved in the experiments using IAD.
6. Applying IAD_{0} to astrophysics: the stability of a Sunlike polytrope
Applying Eqs. (15), (33) and (37) to astrophysics is straightforward by just adding the gravitational acceleration to the momentum equation. We approached gravity using a multipolar expansion of the force (Hernquist & Katz 1989), up to quadrupole contributions. The 3D structure of a Sunlike polytrope was simulated using a sample of N = 10^{5} particles with equal mass. To achieve the equilibrium, we proceeded in three stages. First, the radial coordinate of each particle was set following the 1D density profile of the polytrope with its angular position determined at random. In the second stage, we allowed the particle sample to relax under the action of pressure and gravity forces but the movement of the particles was constrained to keep their distance to the center constant. In a third stage, we allowed the particle sample free to approach the equilibrium configuration. We followed the third stage using both IAD_{0} and STD schemes with the AV coefficients set to α = 1, β = 2. The main goal of this test is to show that IAD_{0} (and probably other matrix methods) can be used to simulate astrophysical scenarios with good results, excellent conservation properties and with low computational penalty.
Fig. 11 Evolution towards stability of a Sunlike star approached by a polytrope. Calculations were carried out using IAD_{0} and the standard, STD, schemes. Labels IAD_{0} (1) and IAD_{0} (2) refer to calculations with the AV parameters set to α = 1,β = 2 and α = 1.5,β = 3 respectively. 

Open with DEXTER 
In Fig. 11 (upper row) we depict the evolution of density for two particles initially located at the center and the surface of the star respectively. As we can see, these particles display the typical oscillatory behavior in both regions. As expected, the amplitude of the oscillations is small close to the center and larger at the surface. The relaxation towards equilibrium is slower for the tensor method but the structure of the star after t = 2 h of relaxation is very similar. With the exception of a tiny region close to the surface, we verified that the gradient of pressure matched gravity along the polytrope in both calculations. In Fig. 11 (lower row) we show the evolution of the kinetic energy and the instantaneous position of the center of mass during the relaxation process. The evolution is similar in both calculations but the level of kinetic energy is always a bit higher when IAD_{0} is used. The energy conservation after two hours of evolution is rather good for both schemes ≃3 × 10^{5}. In the bottom right of Fig. 11 we show the evolution of the moduli of the center of mass position. We see that momentum is not perfectly conserved during the simulation in both simulations, mainly owing to the multipolar approach to the gravitational force. Nevertheless, the tensor scheme provides a better conservation of momentum than the standard SPH calculation.
The higher level of kinetic energy and the slower relaxation rate towards complete mechanical equilibrium shown by IAD_{0} suggest that for disordered systems the numerical noise could be higher for the tensor method. As shown in Fig. 11, the recalculation of the evolution of the polytrope using larger values of the AV coefficients, α = 1.5,β = 3 led to a significant reduction in the level of spurious kinetic energy, while the evolution of the density at the center and surface of the star remained unaltered. This implies that the numerical noise is the main agent responsible for the excess of kinetic energy seen during the relaxation process.
To explore the performance of the algorithm, the tolerance parameter θ, controlling the accuracy of the multipolar expansion, was varied from 0 ≤ θ ≤ 1, where θ = 0 corresponds to direct particle to particle interaction. This benchmarking analysis is summarized in Fig. 12. We see that for standard values of the tolerance parameter, θ ≃ 0.7, widely used in current SPH calculations, the computational overload is around 40%. This overload rapidly decays as θ → 0, as expected. However, the number of operations in the multipolar calculation of gravity scales as ≃Nlog N (Hernquist & Katz 1989). We therefore expect a reduction in the computational overload with increasing numbers of particles. Nevertheless, we note that all these calculations were performed using a serial code. An invaluable property of SPH is that only a single treewalk is needed to find the neighbors (and calculate gravity if needed). Once this information is available the remaining calculations (density, momentum and energy equations, IAD_{0} terms) can be performed with linkedlists that can be directly parallelized. Following that idea, we found that for a 3D simulation of 2 × 10^{5} particles in a standard 4core desktop computer, the gravity calculation (which remains serial) takes around 30 times more wallclock time than the rest of the calculations together (parallelized with OMP using 4 cores). Knowing this, the computational overload due to the IAD_{0} calculation remains subdominant.
Fig. 12 Benchmarking of IAD_{0} versus STD. Parameter θ is the tolerance parameter assumed in the multipolar calculation of gravity. 

Open with DEXTER 
7. Conclusions
We have presented and verified a novel procedure (IAD) to evaluate gradients in the context of SPH simulations. The mechanism for calculating derivatives relies on an integral approach, which ensures that the derivative of a linear function is exactly obtained, even after transforming integrals to summations. The drawback is the greater algebraic complexity because gradients are estimated using the tensor expression given by Eq. (4). Nevertheless, our new scheme is not significantly more complex than the standard one, easy to implement and includes the classical formulation as a particular case. A shortcoming of our new method is that it leads to nonconservative movement equations. Thus, we have developed and tested a restricted interpretation of the method, referred to as IAD_{0}, which sacrifices exact linear interpolation to ensure a perfect conservation of momentum and energy. The analytical considerations and numerical experiments described in Sects. 2 and 4 strongly indicate that IAD_{0} behaves better than the standard method in computing gradients. The modified momentum equation obeying these derivative rules was developed in Sect. 3, resulting in Eq. (29). Since the movement equation was obtained from the EulerLagrange variational principle assuming isentropic evolution, the ensuing equation was fully conservative and explicitly included the ∇h terms. A conservative energy equation compatible with the momentum equation was also developed. To handle with shocks the scheme was completed by including of the standard artificial viscosity formalism resulting in Eqs. (33) and (37), which in addition to the density equation in Eq. (15), summarizes the mathematical formalism linked to IAD_{0}.
The formulation of SPH using matrix methods based on the variational approach (Bonet & Lok 1999) has been used in CFD to successfully simulate a variety of problems (generally in two dimensions) from fluids to the impact and fracture of solid bodies. Nevertheless, none of these schemes are able to simultaneously achieve renormalization, exact momentum and energy conservation, perfect linear interpolation and implicitly include the gradient of the smoothing length in the equations (Oger et al. 2007). Nevertheless, IAD_{0} is probably the optimal formulation because it fulfills almost all the above requirements with a moderate computational overload.
Four tests were performed in two dimensions to verify the performance of the method (see Table 1), two of them in connection to bidimensional subsonic hydrodynamic evolution (KH and RT instabilities) and the other two related to the description of highly supersonic phenomena such as the wall heating shock and the Sedov tests. In all cases the performance of the new scheme was superior, although in the supersonic numerical experiments the improvement was modest. The tests of the growth of the KelvinHelmholtz and RayleighTaylor instabilities clearly showed the power of the IAD_{0} scheme because no growth at all was seen when the standard scheme was used with small initial perturbations, whereas instabilities were able to grow when IAD_{0} was used. This does not of course mean that the standard SPH technique cannot give a satisfactory answer to these problems but simply states that for the same initial conditions the tensor method seems to be more stable and sensitive to small perturbations. Therefore, these results indicate that simulations using the IAD_{0} scheme are less dependent on the initial setting of the particles. As expected of a tensor method, this advantage could probably be reinforced as the dimensionality increases.
In Sect. 5 we devised and discussed a slightly different formulation, which combines exact linear interpolation and good conservative properties, given by Eqs. (52) and (53). The main point here is that the fully IAD method would be only useful for describing linearized systems with smooth density gradients. Nevertheless, a single parameter allows us to easily switch from IAD to IAD_{0} and in this sense the formulation given in Sect. 5 is more general. As sound waves are the prototype of a linearized system we have used the new scheme to simulate the 2D propagation of acoustic waves. The simulations indicate that there are no particular advantages of using IAD instead of either IAD_{0} or standard SPH, because only the symmetry was better preserved during the evolution of the wave trains.
A preliminary application of IAD_{0} to astrophysics was discussed in Sect. 6 in connection to the stability of a Sunlike star described by a polytropic EOS. The main goal was to demonstrate that matrix SPH methods can be used to handle 3D selfgravitating bodies with satisfactory results and to explore the real computational overload when longrange forces need to be computed. For a reasonable accuracy in the calculation of gravity, we have estimated a computational overload <50% for a serial code, with respect to that of standard SPH, when a multipolar expansion is used to calculate the gravitational force. The overload becomes negligible for a precise particletoparticle calculation of gravity scaling as N^{2}/2.
To summarize we could say that simulations carried out using the SPH scheme obtained with the IAD_{0} approach to the gradients always led to improved results with respect to standard SPH. As the new scheme is both fully conservative and more precise in making interpolations, it could be an alternative to the standard technique in handling systems subjected to small perturbations. This conclusion is supported by the results of our numerical study of the growth of KelvinHelmholtz and RayleighTaylor instabilities. In addition, our simulations of supersonic phenomena also improved when the tensor approach was used. The main drawback of this method is that it increases the computational overload, but one has to keep in mind that there is not always a linear relationship between algebraic complexity and computational charge. This could be true for hydrocodes that incorporate timeconsuming physics or when the SPH algorithm built with IAD_{0} allows longer time steps. Even more, using linkedlists a direct parallelization of the method is possible as the calculations are carried out of the treewalk, keeping the computational overload of the IAD_{0} almost negligible compared to more timeconsuming sections of the code.
Although the presented results are encouraging more work needs to be done to confirm and extend the conclusions of our proposal. For the most part, the tests presented in this work to validate the scheme were carried out in 2D boxes using wellordered initial models. The simulation of realistic astrophysical scenarios generally involves, however, a quite different numerical setting. Many of these calculations have to be performed in 3D using a randomlike initial particle distribution and incorporating the longrange gravitational force. A detailed comparative analysis of the ability of IAD_{0} and standard SPH to cope with these scenarios is beyond the scope of the present work and is left to a forthcoming publication. Nonetheless, we can suggest a couple of areas of difficulty in the tensor formulation that will probably come up in astrophysical 3D applications of the method: 1) As suggested in Sect. 6, the tensor method might be more sensitive to the disorder of the particles and display a higher level of numerical noise than the standard scheme. 2) Free boundary conditions could be more difficult to handle using IAD_{0} because the simplification hypothesis assumed in Eq. (14) does not hold at the edge of the system. The difficulty
in simulating sharp boundaries with matrix methods is a wellknown problem of MLS schemes. According to Oger et al. (2007), it can be solved by taking special conditions at the limits of the system. This is not, however, a big concern in astrophysics because astronomical bodies never have abrupt boundaries. Thus, the stability test of the Sunlike star discussed in Sect. 6 does not reveal that the particles located close to the surface have peculiar behaviors. Although the initial imbalance between gravity and pressure gradient is more pronounced than in the standard formalism, the ensuing evolution towards equilibrium is not much different.
Acknowledgments
The authors acknowledge fruitful discussions with Antonio Relaño and Fayyaz Ahmad regarding the SPH technique. This work has been funded by the Spanish MEC grants AYA201015685, AYA200804211C02C01 and the Swiss Platform for HighPerformance and HighProductivity Computing within the supernova project. It was also supported by DURSI of the Generalitat de Catalunya. The rendered SPH plots were made using the freely available SPLASH code (Price 2007).
References
 Abel, T. 2011, MNRAS, 413, 271 [NASA ADS] [CrossRef] (In the text)
 Brookshaw, L. 1985, Proc. Astron. Soc. Aust., 6, 207 [NASA ADS] (In the text)
 Bonet, J., & Lok, T. S.L 1999, Comput. Meth. Appl. Mech. Eng., 180, 97 [NASA ADS] [CrossRef] (In the text)
 Cabezón, R. M., GarcíaSenz, D., & Relaño, A. 2008, J. Comput. Phys., 227, 8523 [NASA ADS] [CrossRef] (In the text)
 Dilts, G. A. 1999, Int. J. Numer. Methods, 44, 1115 [CrossRef] (In the text)
 Gingold, R. A., & Monaghan, J. J., 1977, MNRAS, 181, 375 [NASA ADS] [CrossRef] (In the text)
 Hernquist, L., & Katz, N. 1989, ApJS, 70, 419 [NASA ADS] [CrossRef] (In the text)
 Layzer, D. 1955, ApJ, 122, 1 [NASA ADS] [CrossRef] (In the text)
 Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] (In the text)
 Monaghan, J. J. 1992, ARA&A, 365, 199 (In the text)
 Monaghan, J. J. 2005, Rep. Prog. Phys., 68, 1703 [NASA ADS] [CrossRef] (In the text)
 Morris, J., & Monaghan, J. J. 1997, J. Comput. Phys., 136, 41 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Noh, W. F. 1987, J. Comput. Phys., 72, 78 [NASA ADS] [CrossRef] (In the text)
 Oger, G., Doring, M., Alessandrini, B., & Ferrant, P. 2007, J. Comput. Phys., 225, 1472 [NASA ADS] [CrossRef] (In the text)
 Price, D. 2007, PASA, 24, 159 [NASA ADS] [CrossRef] (In the text)
 Price, D. 2008, J. Comput. Phys., 227, 10040 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Robertson, B. E., Kravtsov, A. V., Gnedin, N. Y., Abel, T., & Rudd, D. H. 2010, MNRAS, 401, 2463 [NASA ADS] [CrossRef] (In the text)
 Rosswog, S. 2009, New Astron. Rev., 53, 78 [NASA ADS] [CrossRef] (In the text)
 Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academic Press Inc.) (In the text)
 Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649 [NASA ADS] [CrossRef] (In the text)
 Youngs, D. 1984, Physica D, 12, 32 [NASA ADS] [CrossRef] (In the text)
All Tables
Conservation properties of the computed models: KelvinHelmholtz (KH), RayleighTaylor (RT), the wall heating shock (WH) and Sedov (SED).
All Figures
Fig. 1 Relative error ϵ in the first derivative of density ρ(x) = (1 + x) g cm^{3} as a function of the smoothing length (in interparticle units) calculated using both tensor and standard SPH schemes. Crosses are for direct derivative estimation. Upper left is for the error when the Gaussian kernel is used. The upper right picture shows the contribution of several error sources to the total error ϵ. Bottom pictures are the same but for the cubicspline kernel. 

Open with DEXTER  
In the text 
Fig. 2 Relative error ϵ in the first derivative of linear functions calculated in 2D as a function of the smoothing length. Upper rows depict the value of ϵ for ρ(x,y) = (1 + x) g cm^{3} obtained using the Gaussian (left) and the cubic spline (right) kernels. Bottom rows are the same but for the function g cm^{3}. Error profiles with crosses are for the direct derivative estimation. 

Open with DEXTER  
In the text 
Fig. 3 Evolution of magnitude of an inhomogeneous 2D system, initially in hydrostatic equilibrium, calculated using the different SPH schemes mentioned in the text. 

Open with DEXTER  
In the text 
Fig. 4 Evolution of the KelvinHelmholtz instability. The snapshots show the color map of density at times t = 0.1,1.0,2.0 and 3.0 s for methods IAD_{0} (first row) and standard SPH (second row), with Δv_{y} = 0.1 cm s^{1}. Times t = 1.0,3.0,4.0 and 5.0 s for methods IAD_{0} (third row) and standard SPH (fourth row), with Δv_{y} = 0.01 cm s^{1}. 

Open with DEXTER  
In the text 
Fig. 5 Evolution of the absolute deviation from initial linear momentum during the development of the KelvinHelmholtz instability. Figure on the left is for the x component of total momentum calculated using IAD_{0} and standard (STD) schemes, whereas figure on the right is the same for the y component. 

Open with DEXTER  
In the text 
Fig. 6 Development of the RayleighTaylor instability for Atwood number 1/3 calculated using the IAD_{0} scheme. The snapshots show the color map of density at times t = 0.4,4.2,5.1 and 5.7 s. 

Open with DEXTER  
In the text 
Fig. 7 Evolution of the center of mass of the spikes during the growth of the RayleighTaylor instability depicted in Fig. 6. The velocity of the center of mass of the bubbles, v_{b}, and spikes, v_{s}, is also shown. For t > 4 s, a limiting speed is reached. 

Open with DEXTER  
In the text 
Fig. 8 Density and velocity profiles during the wall heating shock test. Figure on the left is for IAD_{0} calculation, whereas figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 16 g cm^{3}. 

Open with DEXTER  
In the text 
Fig. 9 Density and velocity profiles during the Sedov test. Figure on the left is for IAD_{0} calculation, whereas the figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 4 g cm^{3}. 

Open with DEXTER  
In the text 
Fig. 10 Acoustic wave profiles at t = 0.35 s calculated using IAD, IAD_{0} and standard SPH (STD). The waves were generated by the periodic displacement of a circular piston of size 0.035 cm located at the center of the lattice. The continuum line is the analytical solution calculated taking a sound speed of c_{s} = 1.29 cm s^{1}. 

Open with DEXTER  
In the text 
Fig. 11 Evolution towards stability of a Sunlike star approached by a polytrope. Calculations were carried out using IAD_{0} and the standard, STD, schemes. Labels IAD_{0} (1) and IAD_{0} (2) refer to calculations with the AV parameters set to α = 1,β = 2 and α = 1.5,β = 3 respectively. 

Open with DEXTER  
In the text 
Fig. 12 Benchmarking of IAD_{0} versus STD. Parameter θ is the tolerance parameter assumed in the multipolar calculation of gravity. 

Open with DEXTER  
In the text 