© ESO, 2012

1. Introduction

The hydrodynamic method known as smoothed particle hydrodynamics (SPH) is a grid-less 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 large-scale 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 well-identified 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), non-standard 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 Euler-Lagrange 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 three-dimensional (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 $\mathit{I}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{V}}\left[\mathit{f}\mathrm{\left(}{{r}}^{\mathbf{\prime }}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\right]\mathrm{(}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{)}\mathit{W}\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathit{,h}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{3}\mathit{,}$(1)where f(r) is any differentiable function and W(|r′ − r|,h) is a spherically symmetric interpolating kernel, usually a sharp-like 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 $\mathit{f}\mathrm{\left(}{{r}}^{\mathbf{\prime }}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{\nabla }{f}\mathrm{·}\mathrm{(}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{)}\mathrm{+}\mathrm{higher}\mathrm{order}\mathrm{terms}$(2)and writing $\mathrm{(}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{)}\mathrm{=}\mathrm{(}{\mathit{x}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{-}{\mathit{x}}_{\mathrm{1}}\mathrm{)}{i}\mathrm{+}\mathrm{(}{\mathit{x}}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{-}{\mathit{x}}_{\mathrm{2}}\mathrm{)}{j}\mathrm{+}\mathrm{(}{\mathit{x}}_{\mathrm{3}}^{\mathrm{\prime }}\mathrm{-}{\mathit{x}}_{\mathrm{3}}\mathrm{)}{k}\mathit{,}$(3)Eq. (1) can be expressed as a matrix equation: $\left[\begin{array}{c}\\ \mathit{\partial f}\mathit{/}\mathit{\partial }{\mathit{x}}_{\mathrm{1}}\\ \mathit{\partial f}\mathit{/}\mathit{\partial }{\mathit{x}}_{\mathrm{2}}\\ \mathit{\partial f}\mathit{/}\mathit{\partial }{\mathit{x}}_{\mathrm{3}}\end{array}\right]\mathrm{=}{\left[\begin{array}{c}\\ {\mathit{\tau }}_{\mathrm{11}}& {\mathit{\tau }}_{\mathrm{12}}& {\mathit{\tau }}_{\mathrm{13}}\\ {\mathit{\tau }}_{\mathrm{21}}& {\mathit{\tau }}_{\mathrm{22}}& {\mathit{\tau }}_{\mathrm{23}}\\ {\mathit{\tau }}_{\mathrm{31}}& {\mathit{\tau }}_{\mathrm{32}}& {\mathit{\tau }}_{\mathrm{33}}\end{array}\right]}^{-1}\left[\begin{array}{c}\\ {\mathit{I}}_{\mathrm{1}}\\ {\mathit{I}}_{\mathrm{2}}\\ {\mathit{I}}_{\mathrm{3}}\end{array}\right]\mathit{,}$(4)where ${\mathit{\tau }}_{\mathit{ij}}\mathrm{=}\mathrm{\int }\mathrm{(}{\mathit{x}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}{\mathit{x}}_{\mathit{i}}\mathrm{)}\mathrm{(}{\mathit{x}}_{\mathit{j}}^{\mathrm{\prime }}\mathrm{-}{\mathit{x}}_{\mathit{j}}\mathrm{)}\mathit{W}\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathit{,h}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{3}\mathrm{;}\mathit{i,j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{3}$(5)and ${\mathit{I}}_{\mathit{k}}\mathrm{=}\mathrm{\int }\left[\mathit{f}\mathrm{\left(}{{r}}^{\mathbf{\prime }}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\right]\mathrm{(}{\mathit{x}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}{\mathit{x}}_{\mathit{k}}\mathrm{)}\mathit{W}\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathit{,h}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{3}\mathrm{;}\mathit{k}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{3.}$(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, τ₁₁ = τ₂₂ = τ₃₃ and τ_ij = τ_ji = 0, ij. Taking for instance the Gaussian kernel ${\mathit{W}}^{\mathit{G}}\mathrm{\left(}\mathit{u,h}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\pi }}^{\mathit{n}\mathit{/}\mathrm{2}}}\exp \mathrm{(}\mathrm{-}{\mathit{u}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(7)where n is the dimensionality of space and $\mathit{u}\mathrm{=}\frac{\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathrm{\right)}}{\mathit{h}}$, it leads to ${\mathit{\tau }}_{\mathit{ii}}\mathrm{=}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{2}}$. Therefore, for the Gaussian kernel the gradient of f reduces to $\mathrm{\nabla }{f}\mathrm{=}\mathrm{\int }\left[\mathit{f}\mathrm{\left(}{{r}}^{\mathbf{\prime }}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\right]\mathrm{\nabla }{W}\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathit{,h}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{3}\mathit{,}$(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 $\mathrm{\nabla }{f}\mathrm{=}\mathrm{\int }\mathit{f}\mathrm{\left(}{{r}}^{\mathbf{\prime }}\mathrm{\right)}\mathrm{\nabla }{W}\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathit{,h}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{3}\mathit{,}$(9)which in SPH parlance becomes ${\mathrm{\nabla }}_{\mathit{a}}{f}\mathrm{=}\sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}{\mathit{f}}_{\mathit{b}}\mathrm{\nabla }{W}\mathit{ab}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathit{,}$(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: ${\left(\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{x}}\right)}_{\mathit{a}}\mathrm{=}\frac{{\sum }_{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathrm{(}{\mathit{f}}_{\mathit{b}}\mathrm{-}{\mathit{f}}_{\mathit{a}}\mathrm{)}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}}{{\sum }_{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}{\mathrm{)}}^{\mathrm{2}}{\mathit{W}}_{\mathit{ab}}}\mathit{,}$(11)taking f_a = px_a + q,f_b = px_b + q the expression above gives ${{}^{\mathrm{\right(}}\frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{x}}^{\mathrm{\left)}}}_{\mathit{a}}\mathrm{=}\mathit{p}$.

In two dimensions, the explicit solution for Eq. (4) is written as $\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{x}}_{\mathrm{1}}}\mathrm{=}{\mathit{c}}_{\mathrm{11}}{\mathit{I}}_{\mathrm{1}}\mathrm{+}{\mathit{c}}_{\mathrm{12}}{\mathit{I}}_{\mathrm{2}}\mathrm{;}\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{x}}_{\mathrm{2}}}\mathrm{=}{\mathit{c}}_{\mathrm{21}}{\mathit{I}}_{\mathrm{1}}\mathrm{+}{\mathit{c}}_{\mathrm{22}}{\mathit{I}}_{\mathrm{2}}\mathit{,}$(12)where $\begin{array}{c}\\ {\mathit{c}}_{\mathrm{11}}\mathrm{=}{\left({\mathit{\tau }}_{\mathrm{11}}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{12}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{22}}}\right)}^{-1}\mathrm{;}& {\mathit{c}}_{\mathrm{12}}\mathrm{=}\mathrm{-}\mathrm{(}{\frac{{\mathit{\tau }}_{\mathrm{12}}}{{\mathit{\tau }}_{\mathrm{22}}}}^{\mathrm{)}}{\left({\mathit{\tau }}_{\mathrm{11}}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{12}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{22}}}\right)}^{-1}\mathrm{;}\\ {\mathit{c}}_{\mathrm{21}}\mathrm{=}\mathrm{-}\mathrm{(}{\frac{{\mathit{\tau }}_{\mathrm{12}}}{{\mathit{\tau }}_{\mathrm{11}}}}^{\mathrm{)}}{\left({\mathit{\tau }}_{\mathrm{22}}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{12}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{11}}}\right)}^{-1}\mathrm{;}& {\mathit{c}}_{\mathrm{22}}\mathrm{=}{\left({\mathit{\tau }}_{\mathrm{22}}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{12}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{11}}}\right)}^{-1}\mathit{.}\end{array}$(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 $\begin{array}{ccc}\mathit{I}\mathrm{\left(}{r}\mathrm{\right)}& \mathrm{=}& \int \left[\mathit{f}\mathrm{\left(}{{r}}^{\mathbf{\prime }}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\right]\mathrm{(}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{)}\mathit{W}\mathrm{\left(}\mathrm{\right|}{{r}}^{\mathbf{\prime }}\mathrm{-}{r}\mathrm{\left|}\mathit{,h}\mathrm{\right)}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathrm{3}\\ & \mathrm{\simeq }& \sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathit{f}\mathrm{\left(}{{r}}_{{b}}\mathrm{\right)}\mathrm{(}{{r}}_{{b}}\mathrm{-}{{r}}_{{a}}\mathrm{)}\mathit{W}\mathrm{\left(}\mathrm{\right|}{{r}}_{{b}}\mathrm{-}{{r}}_{{a}}\mathrm{\left|}\mathit{,}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathit{.}\end{array}$(14)Thereafter, the formulation that results after substituting Eq. (14) into Eq. (4) is labeled IAD₀. 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₀ 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₀ 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 x-axis. The density was calculated using ${\mathit{\rho }}_{\mathit{a}}\mathrm{=}\sum _{\mathit{b}}{\mathit{m}}_{\mathit{b}}{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}\mathrm{\right|}{{r}}_{{b}}\mathbf{-}{{r}}_{{a}}\mathbf{|}\mathit{,}{\mathit{h}}_{\mathit{a}}\mathrm{)}\mathit{.}$(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 $\begin{array}{ccc}\left(\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{x}}\right)\begin{array}{}\mathrm{STD}\\ \mathit{a}\end{array}& \mathrm{=}& \sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}{\mathit{\rho }}_{\mathit{a}}{\mathrm{\nabla }}_{\mathit{a}}{\mathit{W}}_{\mathit{ab}}\mathrm{+}\mathit{p}\sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}\mathrm{)}{\mathrm{\nabla }}_{\mathit{a}}{\mathit{W}}_{\mathit{ab}}\\ & \mathrm{=}& {\mathit{ϵ}}_{\mathrm{1}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{2}}\mathit{p,}\end{array}$(16)where p is the “exact” derivative and ϵ₁ ≃ 0, ϵ₂ ≃ 1. On the other hand, there is only one source of error, ϵ₃ ≃ 0,  if IAD₀ is used, given by $\left(\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{x}}\right)\begin{array}{}\mathrm{IA}{\mathrm{D}}_{\mathrm{0}}\\ \mathit{a}\end{array}\mathrm{=}\frac{{\sum }_{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}{\mathit{\rho }}_{\mathit{a}}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}}{{\sum }_{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}{\mathrm{)}}^{\mathrm{2}}{\mathit{W}}_{\mathit{ab}}}\mathrm{+}\mathit{p}\mathrm{=}{\mathit{ϵ}}_{\mathrm{3}}\mathrm{+}\mathit{p}\mathit{.}$(17)In Fig. 1 (upper and bottom left) we represent the value of the relative error in the derivative $\mathit{ϵ}\mathrm{=}\mathrm{|}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{x}}\mathrm{-}\mathit{p}\mathrm{|}\mathit{/}\mathit{p}$, with respect to the analytical value p as a function of the smoothing length h, normalized to the inter-particle 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₀ or the standard derivative were used, the error increased appreciably. The error is large when the smoothing length is shorter than the inter-particle 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₀ follows almost exactly the error curve for the gradient of density calculated as a simple quotient $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}}$ for adjacent particles. In this sense, IAD₀ 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 $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}}$ 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 cubic-spline kernel.

The contributions to the total error of ϵ₁,ϵ₂ and ϵ₃ in Eqs. (16) and (17) are shown in the rightmost part of Fig. 1. The contribution of ϵ₁ and ϵ₃ are similar but the error ϵ₂, 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 bi-dimensional system was set using a sample of N = 62500 equally spaced particles on a two-dimensional lattice. The mass of the particles along the x-axis was conveniently modified to reproduce a linear density profile, ρ(x,y) = 1 + x. The first derivative of the density, $\frac{\mathit{\partial \rho }}{\mathit{\partial x}}$ 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, $\mathit{ϵ}\mathrm{=}\mathrm{|}\frac{\mathit{\partial \rho }}{\mathit{\partial x}}\mathrm{-}\mathit{p}\mathrm{|}\mathit{/}\mathit{p}$, 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 $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}}$ (crosses in Fig. 2) remained closer to the results obtained using IAD₀.

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 $\mathit{\rho }\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$ g cm-3. Error profiles with crosses are for the direct $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{r}}_{\mathit{b}}\mathrm{-}{\mathit{r}}_{\mathit{a}}}$ derivative estimation.

Finally, we considered a spherically symmetric linear profile ρ(x,y) = 1 + r, where $\mathit{r}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$. 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 inter-particle 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₀ are always smaller than that of the standard scheme. Again the gradient of density estimated as a simply quotient $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{r}}_{\mathit{b}}\mathrm{-}{\mathit{r}}_{\mathit{a}}}$ is more closely reproduced by IAD₀ than by the other schemes.

3. The momentum and energy SPH equations using IAD₀

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 Euler-Lagrange 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 non-dissipative systems. It is shown below how the tensor approach built using the IAD₀ scheme can also be compatible with the variational principle.

The Lagrangian of the system is evaluated as $\mathit{L}\mathrm{=}\sum _{\mathit{b}}{\mathit{m}}_{\mathit{b}}\left(\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathit{b}}^{\mathrm{2}}\mathrm{-}{\mathit{u}}_{\mathit{b}}\mathrm{\left(}{\mathit{\rho }}_{\mathit{b}}\mathit{,}{\mathit{s}}_{\mathit{b}}\mathrm{\right)}\right)\mathit{,}$(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 Euler-Lagrange equations subjected to the constraint ρh³ = const., lead to the equation of movement (Springel & Hernquist 2002)¹${\mathit{m}}_{\mathit{a}}\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{a}}\end{array}\mathrm{=}\mathrm{-}\sum _{\mathit{b}}{\mathit{m}}_{\mathit{b}}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathit{b}}}{\mathrm{\nabla }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{b}}\mathit{,}$(19)where ${\mathrm{\Omega }}_{\mathit{b}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\partial \rho }\mathit{/}\mathit{\partial h}{\sum }_{\mathit{c}}{\mathit{m}}_{\mathit{c}}\mathit{\partial }{\mathit{W}}_{\mathit{bc}}\mathit{/}\mathit{\partial h}\mathrm{)}$ is a term accounting for the gradient of h. The i-component in Eq. (19) can be written ${\mathit{m}}_{\mathit{a}}\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{i,a}}\end{array}\mathrm{=}\mathrm{-}\sum _{\mathit{b}\ne \mathit{a}}{\mathit{m}}_{\mathit{b}}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathit{b}}}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{b}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,a}}}\mathrm{-}{\mathit{m}}_{\mathit{a}}\frac{{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathit{a}}}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{a}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,a}}}\mathrm{·}$(20)An estimation of the density gradients using the tensor Eq. (4) is $\mathrm{\nabla }\mathit{\rho }\mathrm{=}\mathrm{𝒞}{I}\mathit{,}$(21)where . Elements of matrix are those defined in Eq. (5) after changing integrals to summations. For particle a, they are defined to be ${\mathit{\tau }}_{\mathit{ij,a}}\mathrm{=}\sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathrm{(}{\mathit{x}}_{\mathit{i,b}}\mathrm{-}{\mathit{x}}_{\mathit{i,a}}\mathrm{)}\mathrm{(}{\mathit{x}}_{\mathit{j,b}}\mathrm{-}{\mathit{x}}_{\mathit{j,a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathit{,}$(22)and the j-component of vector I in the IAD₀ approach is ${\mathit{I}}_{\mathit{j,a}}\mathrm{=}\sum _{\mathit{b}}{\mathit{m}}_{\mathit{b}}\mathrm{(}{\mathit{x}}_{\mathit{j,b}}\mathrm{-}{\mathit{x}}_{\mathit{j,a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathit{.}$(23)Thus, for particle a the i-component of density gradient is $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{a}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,a}}}\mathrm{=}& {\sum }_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{d}}{\mathit{c}}_{\mathit{ij,a}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}{\mathit{I}}_{\mathit{j,a}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{=}& \\ & & \sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{d}}\sum _{\mathit{b}}^{{\mathit{n}}_{\mathit{b}}}{\mathit{m}}_{\mathit{b}}{\mathit{c}}_{\mathit{ij,a}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{(}{\mathit{x}}_{\mathit{j,b}}\mathrm{-}{\mathit{x}}_{\mathit{j,a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathit{,}\end{array}$(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$\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{b}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,b}}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{d}}\sum _{\mathit{c}}^{{\mathit{n}}_{\mathit{b}}}{\mathit{m}}_{\mathit{c}}{\mathit{c}}_{\mathit{ij,b}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\mathrm{(}{\mathit{x}}_{\mathit{j,c}}\mathrm{-}{\mathit{x}}_{\mathit{j,b}}\mathrm{)}{\mathit{W}}_{\mathit{bc}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\mathit{.}$(25)To calculate $\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{b}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,a}}}$ we use Eq. (15): $\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{b}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,a}}}\mathrm{=}\mathrm{-}{\mathit{m}}_{\mathit{a}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,b}}}\mathrm{·}$(26)The last term in Eq. (26) also appears during the evaluation of ∇_bρ_b, namely $\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{b}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,b}}}\mathrm{=}\sum _{\mathit{c}}{\mathit{m}}_{\mathit{c}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{bc}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,b}}}\mathrm{=}\mathit{...}\mathrm{+}{\mathit{m}}_{\mathit{a}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,b}}}\mathrm{+}\mathit{...}$(27)Comparing Eqs. (27), (26) and (25), we finally get $\frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{b}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i,a}}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{d}}{\mathit{m}}_{\mathit{a}}{\mathit{c}}_{\mathit{ij,b}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\mathrm{(}{\mathit{x}}_{\mathit{j,b}}\mathrm{-}{\mathit{x}}_{\mathit{j,a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\mathit{.}$(28)Substituting Eqs. (24) and (28) into Eq. (20), the i-component of the momentum equation for particle a is given by $\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{i,a}}\end{array}\mathrm{=}\mathrm{-}\sum _{\mathit{b}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{b}}}{\mathit{m}}_{\mathit{b}}\left(\frac{{\mathit{P}}_{\mathit{a}}}{{\mathrm{\Omega }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathrm{\Omega }}_{\mathit{b}}{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\right)\mathit{,}$(29)where $\begin{array}{ccc}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}& & \\ {\mathrm{𝒜}}_{\mathit{i,ab}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}& & \end{array}$It is then straightforward to derive the appropriate energy equation ${\left(\frac{\mathrm{d}\mathit{u}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{a}}\mathrm{=}\frac{{\mathit{P}}_{\mathit{a}}}{{\mathrm{\Omega }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}\sum _{\mathit{b}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{b}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{d}}{\mathit{m}}_{\mathit{b}}\mathrm{(}{\mathit{v}}_{\mathit{i,a}}\mathrm{-}{\mathit{v}}_{\mathit{i,b}}\mathrm{)}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathit{.}$(32)Since the magnitudes , defined by Eqs. (30) and (31), are antisymmetric, the conservation properties of the IAD₀ 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 $\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{i,a}}\end{array}\mathrm{=}\mathrm{-}\sum _{\mathit{b}\mathrm{=}\mathrm{1}}^{\mathit{nb}}{\mathit{m}}_{\mathit{b}}\left(\frac{{\mathit{P}}_{\mathit{a}}}{{\mathrm{\Omega }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathrm{\Omega }}_{\mathit{b}}{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\mathrm{+}{\mathrm{\Pi }}_{\mathit{ab}}\begin{array}{}\mathrm{􏽥}\\ {\mathrm{𝒜}}_{\mathit{i,ab}}\end{array}\right)\mathit{,}$(33)where Π_ab accounts for the viscous pressure ${\mathrm{\Pi }}_{\mathit{ab}}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{-}\mathit{\alpha }{\mathit{c}}_{\mathit{ab}}{\mathit{\mu }}_{\mathit{ab}}\mathrm{+}\mathit{\beta }{\mathit{\mu }}_{\mathit{ab}}^{\mathrm{2}}}{\mathit{\rho ̅}\mathit{ab}}& \mathrm{for}& {r}\mathit{ab}\mathrm{·}{v}\mathit{ab}\mathit{<}\mathrm{0}\mathit{,}\\ \mathrm{0}\mathrm{otherwise}\mathit{,}\end{array}$(34)where the symbols have their usual meaning and r_ab means r_a − r_b. The coefficient μ_ab is ${\mathit{\mu }}_{\mathit{ab}}\mathrm{=}\frac{\mathit{h̅}\mathit{ab}{r}\mathit{ab}\mathrm{·}{v}\mathit{ab}}{{\mathit{r}}_{\mathit{ab}}^{\mathrm{2}}\mathrm{+}\mathrm{0.01}\mathit{h̅}\begin{array}{}\mathrm{2}\\ \mathit{ab}\end{array}}\mathrm{·}$(35)To preserve momentum conservation, we evaluate the arithmetic mean of $\begin{array}{}\mathrm{􏽥}\\ {\mathrm{𝒜}}_{\mathit{i,ab}}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{+}{{\mathrm{𝒜}}^{\mathrm{\prime }}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\right]\mathit{.}$(36)The inclusion of artificial viscosity in the energy equation leads to ${\left(\frac{\mathrm{d}\mathit{u}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{a}}\mathrm{=}\sum _{\mathit{b}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{b}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{d}}{\mathit{m}}_{\mathit{b}}\mathrm{(}{\mathit{v}}_{\mathit{i,a}}\mathrm{-}{\mathit{v}}_{\mathit{i,b}}\mathrm{)}\left(\frac{{\mathit{P}}_{\mathit{a}}}{{\mathrm{\Omega }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{+}\frac{{\mathrm{\Pi }}_{\mathit{ab}}}{\mathrm{2}}\begin{array}{}\mathrm{􏽥}\\ {\mathrm{𝒜}}_{\mathit{i,ab}}\end{array}\right)\mathit{.}$(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 ρ₀ = 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₀ = 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₀ 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₀ 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 ${\mathit{\sigma }}_{\mathit{P}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{=}\sqrt{\frac{{\sum }_{\mathit{b}}\mathrm{(}{\mathit{P}}_{\mathit{b}}\mathrm{-}{\mathit{P}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathit{N}}}$ of an inhomogeneous 2D system, initially in hydrostatic equilibrium, calculated using the different SPH schemes mentioned in the text.

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 tree-walk 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 long-range 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₀ 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₀.

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 Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) hydrodynamic instabilities, whose evolution always remains subsonic and, as we demonstrate, the fully conservative IAD₀ 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 self-adaptive 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₀ formulation as it comes from the Lagrange-Euler 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 point-like 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.

4.1. Kelvin-Helmholtz instability

The Kelvin-Helmholtz instability is a well-known 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 ρ₁ and ρ₂ 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 (ρ₁) moving in a low-density medium (ρ₂) 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 $\mathit{f}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{A}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\exp \frac{\mathrm{2}\mathrm{(}\mathit{y}\mathrm{-}\mathrm{0.25}\mathrm{)}}{\mathrm{\Delta }\mathit{y}}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\exp \frac{\mathrm{2}\mathrm{(}\mathrm{0.75}\mathrm{-}\mathit{y}\mathrm{)}}{\mathrm{\Delta }\mathit{y}}}\mathit{,}$(38)where A is a normalization constant and Δy = 0.05 cm. Therefore, the density profile was given by $\mathit{\rho }\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}\mathrm{-}{\mathit{\rho }}_{\mathrm{2}}\mathrm{)}\mathit{f}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathit{,}$(39)where ρ₁ = 2 g cm^-3 and ρ₂ = 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 ${\mathit{v}}_{\mathit{x}}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{1}}\mathrm{-}{\mathit{v}}_{\mathrm{2}}\mathrm{)}\mathit{f}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathit{,}{\mathit{v}}_{\mathit{y}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{\Delta }{\mathit{v}}_{\mathit{y}}\sin \mathrm{\left(}\mathit{n\pi x}\mathrm{\right)}\mathit{,}$(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₁ = 0.5 cm s^-1 and v₂ = −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 Kelvin-Helmholtz instability at different times for the calculation using IAD₀ (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 well-defined 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₀  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 Kelvin-Helmholtz instability. The snapshots show the color map of density at times t = 0.1,1.0,2.0 and 3.0 s for methods IAD0 (first row) and standard SPH (second row), with Δvy = 0.1 cm s-1. Times t = 1.0,3.0,4.0 and 5.0 s for methods IAD0 (third row) and standard SPH (fourth row), with Δvy = 0.01 cm s-1.

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₀ 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₀ and STD. It is clear that in the standard formulation the instability was unable to grow, while it does using IAD₀. 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 well-known problem, called pairing or tensile-instability, where particles get “stuck” if ${\mathit{r}}_{\mathit{ij}}\mathit{<}\frac{\mathrm{2}}{\mathrm{3}}\mathit{h}$ owing to an unstable stress-strain 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 non-physical 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₀ 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₀.

Fig. 5 Evolution of the absolute deviation from initial linear momentum during the development of the Kelvin-Helmholtz instability. Figure on the left is for the x component of total momentum calculated using IAD0 and standard (STD) schemes, whereas figure on the right is the same for the y component.

4.2. Rayleigh-Taylor instability

The simulation of the growth of the Rayleigh-Taylor 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 non-linear 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 $\mathit{\rho }\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathit{y}\mathrm{\ge }\mathrm{\Delta }\mathit{y}\\ & \mathrm{-}\mathrm{\Delta }\mathit{y}\mathit{<}\mathit{y}\mathit{<}\mathrm{\Delta }\mathit{y}\\ & \mathit{y}\mathrm{\le }\mathrm{-}\mathrm{\Delta }\mathit{y}\end{array}$(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 g cm^-3 and ρ₂ = 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 $\mathit{P}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathit{y}\mathrm{\ge }\mathrm{\Delta }\mathit{y}\\ & \\ & \\ & \mathrm{-}\mathrm{\Delta }\mathit{y}\mathit{<}\mathit{y}\mathit{<}\mathrm{\Delta }\mathit{y}\\ & \mathit{y}\mathrm{\le }\mathrm{-}\mathrm{\Delta }\mathit{y}\end{array}$(42)where g is the gravity, here of value g = 0.5 cm s^-2. With this setting, the speed of sound at the inter-phase 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) ${\mathit{v}}_{\mathit{y}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Delta }{\mathit{v}}_{\mathit{y}}}{\mathrm{4}}\left\{\mathrm{1}\mathrm{+}\cos \left[\mathrm{8}\mathit{\pi }\left(\mathit{x}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}}\right)\right]\right\}\left\{\mathrm{1}\mathrm{+}\cos \left[\mathrm{5}\mathit{\pi }\left(\mathit{y}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\right)\right]\right\}\mathit{,}$(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₀. 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 Kelvin-Helmholtz 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 Rayleigh-Taylor instability for Atwood number 1/3 calculated using the IAD0 scheme. The snapshots show the color map of density at times t = 0.4,4.2,5.1 and 5.7 s.

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) = ξ₀exp[Γ(t − t₀)], roughly until the initial perturbation ξ₀ at t₀ (here we took t₀ = 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 ${\mathrm{\Gamma }}_{\mathrm{t}}\mathrm{=}\sqrt{{\mathit{A}}_{\mathrm{t}}\mathit{kg}}$, where A_t = (ρ₂ − ρ₁)/(ρ₂ + ρ₁) 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² 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), ${\mathit{v}}_{\mathit{b}}\mathrm{=}\mathrm{0.51}\sqrt{{\mathit{A}}_{\mathrm{t}}\mathit{gl}}$ 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 Rayleigh-Taylor instability depicted in Fig. 6. The velocity of the center of mass of the bubbles, vb, and spikes, vs, is also shown. For t > 4 s, a limiting speed is reached.

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 well-known 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 non-physical 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₀ 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 well-developed. 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 post-shock 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 IAD0 calculation, whereas figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 16 g cm-3.

4.4. Sedov test

In the Sedov test, the evolution of a shock wave front born as a consequence of a point-like 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 self-similar 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 smoothing-length. 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 post-shock 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 IAD0 calculation, whereas the figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 4 g cm-3.

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 $\mathit{P}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{P}}_{\mathrm{1}}\mathrm{-}{\mathit{P}}_{\mathrm{2}}\mathrm{)}\exp \left[\frac{\mathrm{-}{\mathit{r}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right]\mathit{,}$(44)where P₁ and P₂ are the pressures in zones left and right of the step and σ sets the width of the pressure decay. In the present simulations, P₁ = 10⁴ dyn cm^-2 and P₂ = 1 dyn cm^-2 and σ² = 36 cm², 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 self-similar 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 post-shock zone is smoother in the calculation using IAD₀. The velocity profile is also similar but the tensor calculation showed fewer oscillations and a more ordered behavior in the post-shock 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₀ 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 $\begin{array}{}\mathbf{¨}\\ {r}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathit{P}\mathrm{=}\mathrm{-}\left[\mathrm{\nabla }\left(\frac{\mathit{P}}{\mathit{\rho }}\right)\mathrm{+}\frac{\mathit{P}}{{\mathit{\rho }}^{\mathrm{2}}}\mathrm{\nabla }\mathit{\rho }\right]\mathit{,}$(45)which, using Eq. (10), can be easily accommodated within the SPH formalism $\begin{array}{}\mathbf{¨}\\ {r}\end{array}\mathit{a}\mathrm{=}\mathrm{-}\sum _{\mathit{b}}{\mathit{m}}_{\mathit{b}}\left(\frac{{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}}\right)\mathrm{\nabla }\begin{array}{}\mathrm{˜}\\ {\mathit{W}}_{\mathit{ab}}\end{array}\mathit{.}$(46)This equation can be deduced independently using the variational principle and, provided the kernel is symmetrized $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{W}}_{\mathit{ab}}\end{array}\mathrm{=}\mathrm{0.5}\mathrm{\left(}{\mathit{W}}_{\mathit{ab}}\mathrm{\right(}{\mathit{h}}_{\mathit{a}}\mathrm{)}\mathrm{+}{\mathit{W}}_{\mathit{ab}}\mathrm{(}{\mathit{h}}_{\mathit{b}}\mathrm{\left)}\mathrm{\right)}$, it conserves linear and angular momentum exactly. If instead of Eq. (10) we use Eq. (8) to compute Eq. (45), the resulting equation is $\begin{array}{ccc}\begin{array}{}\mathbf{¨}\\ {r}\end{array}\mathit{a}\mathrm{=}& & \mathrm{-}\left[\mathrm{\nabla }\left(\frac{\mathit{P}}{\mathit{\rho }}\right)\mathrm{+}\frac{\mathit{P}}{{\mathit{\rho }}^{\mathrm{2}}}\mathrm{\nabla }\mathit{\rho }\right]\\ & & \mathrm{=}\mathrm{-}\sum _{\mathit{b}}{\mathit{m}}_{\mathit{b}}\left(\frac{{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}}\right)\mathrm{\nabla }\begin{array}{}\mathrm{˜}\\ {\mathit{W}}_{\mathit{ab}}\end{array}\mathrm{+}\mathrm{2}{\mathit{P}}_{\mathit{a}}\sum _{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{b}}}\mathrm{\nabla }\begin{array}{}\mathrm{˜}\\ {\mathit{W}}_{\mathit{ab}}\end{array}\mathit{.}\end{array}$(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 / ρ² in Eq. (47) assuming a smooth density gradient $\frac{{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}\mathrm{\simeq }\frac{{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{b}}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{\Delta }\mathit{\rho }}{{\mathit{\rho }}_{\mathit{b}}}\right)\mathit{,}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}}\mathrm{\simeq }\frac{{\mathit{P}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{b}}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{\Delta }\mathit{\rho }}{{\mathit{\rho }}_{\mathit{a}}}\right)\mathit{,}$(48)where Δρ = ρ_a − ρ_b and |Δρ| / ρ ≪ 1. Inserting Eq. (48) into Eq. (47) and neglecting higher order terms, the acceleration is given by $\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{a}}\end{array}\mathrm{=}\mathrm{-}\frac{{\sum }_{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\left(\frac{{\mathit{P}}_{\mathit{b}}\mathrm{-}{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}}\right)\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}\mathrm{)}{\mathit{W}}_{\mathit{ab}}}{{\sum }_{\mathit{b}}\frac{{\mathit{m}}_{\mathit{b}}}{{\mathit{\rho }}_{\mathit{b}}}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}{\mathrm{)}}^{\mathrm{2}}{\mathit{W}}_{\mathit{ab}}}\mathrm{·}$(49)If we assume a linear pressure profile within the kernel domain of the particle ${\mathit{P}}_{\mathit{b}}\mathrm{=}{\mathit{P}}_{\mathit{a}}\mathrm{+}{\left(\frac{\mathit{\partial P}}{\mathit{\partial x}}\right)}_{\mathit{a}}\mathrm{(}{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}\mathrm{)}\mathit{,}$(50)the acceleration becomes $\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{a}}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathit{a}}}{\left(\frac{\mathit{\partial P}}{\mathit{\partial x}}\right)}_{\mathit{a}}\mathit{,}$(51)which is the one-dimensional Newton equation, thus leading to an exact value for the acceleration. In consequence, global momentum is well-preserved 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₀ scheme with a corrective term $\begin{array}{ccc}\begin{array}{}\mathrm{¨}\\ {\mathit{x}}_{\mathit{i,a}}\end{array}\mathrm{=}\mathrm{-}\sum _{\mathit{b}\mathrm{=}\mathrm{1}}^{\mathit{nb}}{\mathit{m}}_{\mathit{b}}& & (\frac{{\mathit{P}}_{\mathit{a}}}{{\mathrm{\Omega }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{P}}_{\mathit{b}}}{{\mathrm{\Omega }}_{\mathit{b}}{\mathit{\rho }}_{\mathit{b}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\\ & & \mathrm{+}{\mathrm{\Pi }}_{\mathit{ab}}\begin{array}{}\mathrm{􏽥}\\ {\mathrm{𝒜}}_{\mathit{i,ab}}\end{array}\mathrm{-}\frac{\mathrm{\Psi }{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{b}}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)})\mathit{,}\end{array}$(52)where 0 ≤ Ψ ≤ 2 is a parameter controlling the strength of the applied correction. For Ψ = 0, the already checked IAD₀ 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) $\begin{array}{ccc}{\left(\frac{\mathrm{d}\mathit{u}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{a}}\mathrm{=}\sum _{\mathit{b}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{b}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{d}}& & {\mathit{m}}_{\mathit{b}}\mathrm{(}{\mathit{v}}_{\mathit{i,a}}\mathrm{-}{\mathit{v}}_{\mathit{i,b}}\mathrm{)}(\frac{{\mathit{P}}_{\mathit{a}}}{{\mathrm{\Omega }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)}\\ & & \mathrm{+}\frac{{\mathrm{\Pi }}_{\mathit{ab}}}{\mathrm{2}}\begin{array}{}\mathrm{􏽥}\\ {\mathrm{𝒜}}_{\mathit{i,ab}}\end{array}\mathrm{\left(}{\mathit{h}}_{\mathit{b}}\mathrm{\right)}\mathrm{-}\frac{\mathrm{\Psi }{\mathit{P}}_{\mathit{a}}}{{\mathit{\rho }}_{\mathit{a}}{\mathit{\rho }}_{\mathit{b}}}{\mathit{v}}_{\mathit{i,a}}{\mathrm{𝒜}}_{\mathit{i,ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{a}}\mathrm{\right)})\mathit{.}\end{array}$(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₀ 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 ρ₀ = 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 ${\mathit{c}}_{\mathrm{s}}\mathrm{=}\sqrt{\mathit{\gamma P}\mathit{/}\mathit{\rho }}\mathrm{=}\mathrm{1.29}$ 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 ${\mathit{v}}_{\mathrm{r}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{r}}^{\mathrm{m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\cos \mathrm{\left(}\mathit{\omega t}\mathrm{\right)}\mathit{,}$(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 ${\mathit{v}}_{\mathrm{r}}^{\mathrm{m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{r}}^{\mathrm{m}}\mathrm{\left(}{\mathit{R}}_{\mathrm{p}}\mathrm{\right)}\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{p}}}\right)\mathrm{·}$(55)The resolution places limits on both the period and the maximum velocity at the piston head ${\mathit{v}}_{\mathrm{r}}^{\mathrm{m}}\mathrm{\left(}{\mathit{R}}_{\mathrm{p}}\mathrm{\right)}$. 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 ${\mathit{v}}_{\mathrm{r}}^{\mathrm{m}}\mathrm{\left(}{\mathit{R}}_{\mathrm{p}}\mathrm{\right)}\mathrm{=}\mathrm{0.16}$ cm s^-1.

Fig. 10 Acoustic wave profiles at t = 0.35 s calculated using IAD, IAD0 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 cs = 1.29 cm s-1.

The results of the simulations for the three SPH approaches: IAD, IAD₀ 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₀ 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 $\mathit{A}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{=}{\mathit{A}}_{\mathit{m}}\mathit{J}\mathrm{\left(}\mathit{u}\mathrm{\right)}\cos \mathrm{\left(}\mathit{wt}\mathrm{\right)}\mathit{,}$(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 ${\mathit{A}}_{\mathrm{m}}\mathrm{=}{\mathit{v}}_{\mathrm{r}}^{\mathrm{m}}\mathrm{\left(}{\mathit{R}}_{\mathrm{p}}\mathrm{\right)}$ 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₀ scheme (${\mathit{c}}_{\mathrm{s}}^{\mathrm{ia}{\mathrm{d}}_{\mathrm{0}}}\mathrm{=}\mathrm{1.285}$ cm s^-1) followed by the standard (${\mathit{c}}_{\mathrm{s}}^{\mathrm{std}}\mathrm{=}\mathrm{1.267}$ cm s^-1), whereas full IAD falls a bit short (${\mathit{c}}_{\mathrm{s}}^{\mathrm{iad}}\mathrm{=}\mathrm{1.25}$ 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₀ or even by the standard scheme. Only the spherical symmetry was more accurately preserved in the experiments using IAD.

6. Applying IAD₀ to astrophysics: the stability of a Sun-like 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 Sun-like polytrope was simulated using a sample of N = 10⁵ 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₀ and STD schemes with the AV coefficients set to α = 1, β = 2. The main goal of this test is to show that IAD₀ (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 Sun-like star approached by a polytrope. Calculations were carried out using IAD0 and the standard, STD, schemes. Labels IAD0 (1) and IAD0 (2) refer to calculations with the AV parameters set to α = 1,β = 2 and α = 1.5,β = 3 respectively.

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₀ 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₀ 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 tree-walk 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₀ terms) can be performed with linked-lists that can be directly parallelized. Following that idea, we found that for a 3D simulation of 2 × 10⁵ particles in a standard 4-core desktop computer, the gravity calculation (which remains serial) takes around 30 times more wall-clock time than the rest of the calculations together (parallelized with OMP using 4 cores). Knowing this, the computational overload due to the IAD₀ calculation remains sub-dominant.

Fig. 12 Benchmarking of IAD0 versus STD. Parameter θ is the tolerance parameter assumed in the multipolar calculation of gravity.

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 non-conservative movement equations. Thus, we have developed and tested a restricted interpretation of the method, referred to as IAD₀, 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₀ 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 Euler-Lagrange 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₀.

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 re-normalization, 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₀ 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 bi-dimensional 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 Kelvin-Helmholtz and Rayleigh-Taylor instabilities clearly showed the power of the IAD₀ 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₀ 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₀ 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₀ 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₀ or standard SPH, because only the symmetry was better preserved during the evolution of the wave trains.

A preliminary application of IAD₀ to astrophysics was discussed in Sect. 6 in connection to the stability of a Sun-like star described by a polytropic EOS. The main goal was to demonstrate that matrix SPH methods can be used to handle 3D self-gravitating bodies with satisfactory results and to explore the real computational overload when long-range 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 particle-to-particle calculation of gravity scaling as N²/2.

To summarize we could say that simulations carried out using the SPH scheme obtained with the IAD₀ 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 Kelvin-Helmholtz and Rayleigh-Taylor 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 time-consuming physics or when the SPH algorithm built with IAD₀ allows longer time steps. Even more, using linked-lists a direct parallelization of the method is possible as the calculations are carried out of the tree-walk, keeping the computational overload of the IAD₀ almost negligible compared to more time-consuming 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 well-ordered 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 random-like initial particle distribution and incorporating the long-range gravitational force. A detailed comparative analysis of the ability of IAD₀ 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₀ 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 well-known 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 Sun-like 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 AYA2010-15685, AYA2008-04211-C02-C01 and the Swiss Platform for High-Performance and High-Productivity 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

All Tables

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 $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{x}}_{\mathit{b}}\mathrm{-}{\mathit{x}}_{\mathit{a}}}$ 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 cubic-spline kernel.

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 $\mathit{\rho }\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$ g cm-3. Error profiles with crosses are for the direct $\frac{{\mathit{\rho }}_{\mathit{b}}\mathrm{-}{\mathit{\rho }}_{\mathit{a}}}{{\mathit{r}}_{\mathit{b}}\mathrm{-}{\mathit{r}}_{\mathit{a}}}$ derivative estimation.

In the text

	Fig. 3 Evolution of magnitude ${\mathit{\sigma }}_{\mathit{P}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{=}\sqrt{\frac{{\sum }_{\mathit{b}}\mathrm{(}{\mathit{P}}_{\mathit{b}}\mathrm{-}{\mathit{P}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathit{N}}}$ of an inhomogeneous 2D system, initially in hydrostatic equilibrium, calculated using the different SPH schemes mentioned in the text.
In the text

	Fig. 4 Evolution of the Kelvin-Helmholtz instability. The snapshots show the color map of density at times t = 0.1,1.0,2.0 and 3.0 s for methods IAD0 (first row) and standard SPH (second row), with Δvy = 0.1 cm s-1. Times t = 1.0,3.0,4.0 and 5.0 s for methods IAD0 (third row) and standard SPH (fourth row), with Δvy = 0.01 cm s-1.
In the text

	Fig. 5 Evolution of the absolute deviation from initial linear momentum during the development of the Kelvin-Helmholtz instability. Figure on the left is for the x component of total momentum calculated using IAD0 and standard (STD) schemes, whereas figure on the right is the same for the y component.
In the text

	Fig. 6 Development of the Rayleigh-Taylor instability for Atwood number 1/3 calculated using the IAD0 scheme. The snapshots show the color map of density at times t = 0.4,4.2,5.1 and 5.7 s.
In the text

	Fig. 7 Evolution of the center of mass of the spikes during the growth of the Rayleigh-Taylor instability depicted in Fig. 6. The velocity of the center of mass of the bubbles, vb, and spikes, vs, is also shown. For t > 4 s, a limiting speed is reached.
In the text

	Fig. 8 Density and velocity profiles during the wall heating shock test. Figure on the left is for IAD0 calculation, whereas figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 16 g cm-3.
In the text

	Fig. 9 Density and velocity profiles during the Sedov test. Figure on the left is for IAD0 calculation, whereas the figure on the right is for the simulation using the standard SPH scheme. Density is normalized to 4 g cm-3.
In the text

	Fig. 10 Acoustic wave profiles at t = 0.35 s calculated using IAD, IAD0 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 cs = 1.29 cm s-1.
In the text

	Fig. 11 Evolution towards stability of a Sun-like star approached by a polytrope. Calculations were carried out using IAD0 and the standard, STD, schemes. Labels IAD0 (1) and IAD0 (2) refer to calculations with the AV parameters set to α = 1,β = 2 and α = 1.5,β = 3 respectively.
In the text

	Fig. 12 Benchmarking of IAD0 versus STD. Parameter θ is the tolerance parameter assumed in the multipolar calculation of gravity.
In the text