Issue 
A&A
Volume 570, October 2014



Article Number  A14  
Number of page(s)  15  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201424260  
Published online  07 October 2014 
Equalizing resolution in smoothedparticle hydrodynamics calculations using selfadaptive sinc kernels
^{1}
Dept. de Física i Enginyeria Nuclear, Universitat Politècnica de
Catalunya,
Compte d’Urgell 187,
08036
Barcelona,
Spain
email:
domingo.garcia@upc.edu
^{2}
Institut d’Estudis Espacials de Catalunya,
Gran Capità 24,
08034
Barcelona,
Spain
^{3}
Departement Physik, Universität Basel,
Klingelbergstrasse 82,
4056
Basel,
Switzerland
Received:
23
May
2014
Accepted:
18
August
2014
Context. The smoothedparticle hydrodynamics (SPH) technique is a numerical method for solving gasdynamical problems. It has been applied to simulate the evolution of a wide variety of astrophysical systems. The method has a secondorder accuracy, with a resolution that is usually much higher in the compressed regions than in the diluted zones of the fluid.
Aims. We propose and check a method to balance and equalize the resolution of SPH between high and lowdensity regions. This method relies on the versatility of a family of interpolators called sinc kernels, which allows increasing the interpolation quality by varying only a single parameter (the exponent of the sinc function).
Methods. The proposed method was checked and validated through a number of numerical tests, from standard onedimensional Riemann problems in shock tubes, to multidimensional simulations of explosions, hydrodynamic instabilities, and the collapse of a Sunlike polytrope.
Results. The analysis of the hydrodynamical simulations suggests that the scheme devised to equalize the accuracy improves the treatment of the postshock regions and, in general, of the rarefacted zones of fluids while causing no harm to the growth of hydrodynamic instabilities. The method is robust and easy to implement with a low computational overload. It conserves mass, energy, and momentum and reduces to the standard SPH scheme in regions of the fluid that have smooth density gradients.
Key words: methods: numerical / hydrodynamics
© ESO, 2014
1. Introduction
The hydrodynamical method known as smoothedparticle hydrodynamics (SPH) is a gridless Lagrangian approach to continuum mechanics devised by Gingold & Monaghan (1977) and Lucy (1977). A key ingredient of the SPH technique is the nature of an interpolating function called the kernel, which is used to estimate the value of different physical magnitudes. Because of the nature of the SPH interpolations, the gradient of any magnitude can be calculated by directly taking the gradient of the kernel, which is an analytically differentiable function. This provides an easy and efficient way of obtaining gradients. In that way it is easy to write the Euler equations of fluid mechanics in terms of the kernel and its derivatives (Monaghan 1992, 2005). Despite its success, SPH still has several weak points, which have recently caused a number of improvements of the technique (Rosswog 2014; Saitoh & Makino 2013; GarcíaSenz et al. 2012; Cabezón et al. 2012; Valdarnini 2012; Dehnen & Aly 2012; Springel 2010a). One of the shortcomings of SPH is that the accuracy in the density evaluation is different in all fluid regions. The resolution in lowdensity regions is typically poorer than in the highdensity regions. In astrophysics, the regions close to the surface of selfgravitating bodies usually have a lower resolution than the interior. Another source of inaccuracy are fluid discontinuities, such as shock waves or sharp boundaries. In this case, the difficulty comes from the inefficacy of the interpolations to keep track of phenomena with a lengthscale lower than the characteristic smoothing length h.
Fig. 1 Left: profile of several sinc (S_{n}, continuum lines) and Wendland kernels (ψ_{l:k}) (rescaled to a common range [0, 1]). Points + (in red), × (in green) and ∗ (in blue) are for ψ_{3:1}, ψ_{4:2} and ψ_{5:3} respectively. Right: same as before, but for sinc and cubic, quartic, and quintic spline kernels. 
It is well known that the standard formulation of SPH is secondorder accurate in h. An interpolating function W = (K/h) f(v) is defined so that the averaged value of density at any point of the system is estimated (in 1D) as (1)where K is a normalization constant and v =  x − x_{0}  /h, where h is a scaling parameter called the smoothinglength. The magnitude ⟨ ρ(x_{0}) ⟩ is the SPH estimation of density at the fluid coordinate x_{0}, and ρ_{0} is the value of the density at that point. The smoothing length h is usually taken as the local resolution. Therefore, the lead dependence of the error in evaluating the density is proportional to h^{2}. In most applications the value of the smoothinglength is directly linked to the density via mass conservation, h ∝ ρ^{− 1 /d}, where d is the dimension of the space, resulting in a selfconsistent Lagrangian description of the dynamics (Springel & Hernquist 2002). However, setting the smoothing length in that way is detrimental to the lowdensity regions, where h becomes large. On another note, the error term in Eq. (1) has a dependence on the second derivative of the smoothed function itself, strictly vanishing only for linear functions. As a consequence, one can expect relevant differences in accuracy along the fluid wherever the second derivative of the function does not vanish.
To this extent, we can use the shape of the kernel as a parameter to gain additional control over the error in Eq. (1) so that the interpolations in density are made with approximately the same relative accuracy in any point of the system. An easy way of achieving this is using the sinc kernels.
In this paper we introduce a novel method to balance the resolution in all regions of the fluid. We take advantage of the features of the oneparametric family of kernels introduced by Cabezón et al. (2008) to set an additional variable (the exponent n of the sinc function, given by Eq. (2)) which, jointly with h, controls the local resolution. We found a mechanism to increase the accuracy of interpolations by changing n when a strong density gradient is detected. In fluid regions with low or moderate gradients, the scheme reduces to the standard (Monaghan 2005), where the resolution is basically set by the smoothing length.
The paper is organized as follows: in Sect. 2, the main features of the proposed family of kernels are reviewed and compared with other existing interpolators. The numerical scheme used to equalize the resolution is explained in Sect. 3, where we also provide some insights into the ability of the algorithm to handle sharp onedimensional density profiles. The basic Euler equations, written in the SPH formalism and incorporating the adaptive sinc kernels, are described in Sect. 4 and in Appendix A. These sections also include the treatment of gradh and gradn terms, as well as details of the practical implementation of the algorithm. In Sect. 5, we check the hydrocode with a variety of standard tests in one, two, and three dimensions. These tests intend to cover different physical processes, such as shocks, instabilities, or selfgravitating bodies. Finally, we outline the main conclusions of our work and comments on the shortcomings of the developed scheme as well as on future lines of improvement in Sect. 6.
Fig. 2 Fourier transform,  ℱ_{3}  of sinc kernels with n = 3,4 and 10 as well as the cubic spline, where the dashed lines indicate the negative portions of the curves. They can be compared with the Fourier transform of several Wendland and spline kernels given in Fig. 2 by Dehnen & Aly (2012). 
Fig. 3 Color map of the logarithm of the error in density estimation in 2D (left) and 3D (right) as a function of the kernel index and number of neighbors (see Sect. 3 for a complete explanation of this diagram). The dashed line in red is the rough critical limit separating the region susceptible to particle pairing (above the line). The blue dotted line denotes the region where the approach of integrals by summations in density calculation is too sensitive to particle distribution (below the line). 
2. Main features of the sinc family of kernels
The sinc family of compact supported kernels was first introduced by Cabezón et al. (2008, 2012) (who called them harmonic kernels) as a way to make the SPH technique more flexible. The ultimate goal was to gather in a unique function the more relevant features of some of the most often used interpolators. The sinc family of kernels is defined as (2)where S_{n}(.) = sinc^{n}(.), and where n is the index of the kernel and B_{n} a normalization constant. The function is a widely known function used in signal analysis and spectral theory. The main advantage of this kernel is that it is able to mimic the behavior of some of the most popular kernels, such as the cubic and quintic splines (Cabezón et al. 2008). By a careful choice of the index n, it can even approach several of the socalled Wendland kernels (Wendland 1995), recently discussed in Dehnen & Aly (2012), which behave optimally in avoiding pairing instability. In Fig. 1 we present the profiles of several sinc, S_{n}, and Wendland ψ_{l:k} kernels with the independent variable rescaled to the range [ 0,1 ]. The profiles of the high order ψ_{3:1}, ψ_{4:2}, and ψ_{5:3} Wendland interpolators can be approached by S_{3.97}, S_{5.12}, and S_{6.41} respectively, whereas several members of the spline family can be approached by S_{3.07} (cubic spline), S_{3.97} (quartic spline), and S_{4.89} (quintic spline). According to Dehnen & Aly (2012), a necessary condition to escape the pairing instability is that the Fourier transform, , (in 3D, where κ is the wavenumber) of the kernel is always definite positive. In Fig. 2 we show ℱ_{3}(S_{n}) for kernel indexes n = 3,6,10. For n = 3 (similar to the cubic spline) the Fourier transform becomes negative at relatively low wavelengths, limiting the number of neighbors in the summations in practical applications. Nevertheless, for n = 6 the dynamical range has been considerably extended and is now similar to that of the quintic spline. For n = 10 the Fourier transform becomes marginally negative at long wavelengths, thus showing a good endurance in front of the pairing.
The relationship between the index of the kernel, n, and the maximum number of neighbors (to elude particle clustering), , can be made more explicit by applying the empirical rule described in Price (2012). Many computer simulations have shown that particle clustering is avoided when the normalized interparticle separation in ordered lattices becomes larger than a critical value, r_{ij}/h ≥ η, with η ≃ 1.2 for the cubic spline (Price 2012). A value η = 1.2 corresponds to neighbors in 2D and 3D, respectively, so that for there is an increasing chance for particle clustering. The second derivative of the sinc,n = 3 (cubicsplinelike kernel), vanishes at η′ = r_{ij}/h = 1.5, which is a factor f = 5 / 4 larger than η = 1.2. To find for any n, we first estimated the point η′ where the second derivative of S_{n} vanishes as a function of the kernel exponent n and calculated the value of η(n) = η′f^{1}. That value of η(n) was then mapped to (i.e., the maximum amount of neighbors that is still resistant to particle clustering). In Fig. 3 we depict the profile of the maximum number of neighbors as a function of the kernel index n in two and three dimensions. For neighbors in 3D it is advisable to take n ≥ 3,4, and 7, respectively, in Eq. (2) to avoid particle pairing.
The minimum number of neighbors is also constrained by the approach of SPH integrals as finite summations. The precise value of is difficult to estimate because it ultimately depends on the specific physical problem. A simple, albeit qualitative, way to set is to numerically calculate the value of the density in a regular bcc lattice as a function of the kernel index n and the number of neighbors, and measure of the relative deviation of the density from the theoretical value. To do this, we built 2D and 3D regular square lattices of unit square and placed N particles at the nodes of the grid. By giving a mass m = 1 /N to each particle, this should result in a uniform density ρ = 1. As a practical criterion for selecting , we calculated the density, ρ_{SPH}, for many pairs as well as the value of σ = (  ρ_{SPH} − 1 ). These combinations with σ ≥ 10^{3} were then considered not a good enough realization of density, and the critical limiting values with σ ≃ 10^{3} were stored. The results of this study is summarized in the dotted blue lines shown in Fig. 3. The most striking feature of these lines is that they have a minimum for kernel indexes n = 4 − 6, which indicates that these sinc interpolators are the best choice to carry out SPH calculations. These results support the current feeling that using the quarter or the quintic spline enhances the convergence of interpolations. It is worth noting that when the pairing and the convergence criteria are combined, an optimal domain in the plane appears, which, for a given kernel index n, restricts the number of neighbors to . It is equally remarkable that the pairing and the convergence lines intersect at n ≃ 3 in 2D and 3D, meaning that a large body of SPH calculations carried out so far with the cubicspline kernel have probably been bordering on an unpredictable and dangerous zone. It should be kept in mind, however, that the path of these lines in the diagram is merely qualitative, and our advice is not to proceed too far into the forbidden regions in practical applications.
The implementation of the family of kernels adds more flexibility to SPH because one can, for example, take a different index n to handle the artificial viscosity terms in the momentum and energy equations or in the heat conduction equation, without changing the number of neighbors of the particle. Another virtue of the family is that by a careful choice of kernel exponents it allows equalizing the interpolation accuracy along the system; this last point is the subject of the present study.
To increase the computational speed, it is recommended to store the value of and its derivative as a function of v, (0 ≤ v ≤ 2) in a table and use a linear Taylor expansion to calculate the value of sinc and other variables of interest (see Appendix A). This allows a fast computation of Eq. (2) and its derivative after the index n is chosen. The value of the normalization constant B_{n} for 2 ≤ n ≤ 12 can be obtained from the following fitting functions: (3)where the values of coefficients b_{0},b_{1},b_{2},b_{3} as a function of the dimensionality are provided in Table 1. These fitting functions are fast to compute and precise up to the fifth significant figure.
3. Using the sinc kernels to enhance interpolations
According to Eq. (1), the leading error ℰ in estimating den sity is (4)where n and h are the kernel index and the smoothing length. In the standard Lagrangian formulation of SPH (Monaghan 2005) the smoothing length is constrained so that the neighboring mass of a given particle always remains constant. This fact leaves the kernel index n as the only free parameter to control the error size given by Eq. (4). In general, increasing the value of n decreases the error, but one cannot increase n arbitrarily without generating too much numerical noise.
It is easy to estimate the relative accuracy achieved with several exponents n using Eq. (4). First we write the error term as (5)with (6)with A = 2,4, and p = 3,4 in 2D and 3D, respectively, and B_{n} is the normalization constant given by Eq. (3). Setting arbitrarily for n = 3, neighbors in 3D and for n = 3, in 2D (i.e., we choose the normalization for the error), and taking into account that , we write (7)For a given kernel index n the integral ℐ_{n} is calculated numerically and the value of the magnitude 2ℰ /  ∇^{2}ρ  as a function of pairs is shown in Fig. 3. Although 2ℰ /  ∇^{2}ρ  is not directly the interpolation error, it does give an estimate of it. Therefore, for the sake of simplicity, we refer to this magnitude as the error in the following. From this figure we see, for example, that a similar accuracy is achieved by different pairs . For example, the S_{3}sinc kernel using neighbors in 3D has a similar leading error term as S_{5} and or S_{10} and . On the whole, computations are fast for low but the results are more sensitive to numerical noise, while the opposite is true for high . A conservative option is to work with a moderate number of neighbors and variable kernel indexes, leaving high kernel exponents to handle only regions with steep gradients. As a default, we took n = 5 and neighbors in the 2D and 3D numerical experiments described in Sect. 5, although other combinations are also feasible. A similar number of neighbors and a high Bspline kernel, quartic (≃ S_{4}) or quintic (≃ S_{5}), was suggested by Valdarnini (2012) as an optimal choice to improve the convergence of hydrodynamic simulations.
Fig. 4 Fitting of several 1D density profiles using the sinc family of kernels with constant and selfadaptive indexes. The upper left panel shows a Gaussian, mountainlike, profile. The exact (e) analytical value is shown in red, the result with the variable (v) kernel index, constant n = 5 and n = 3 in green, blue, and pink, respectively. The light blue line shows the profile of the smoothinglength normalized to its highest value (achieved just at the limits of the system, x = 0 and x = 1). The orange line shows the profile of the kernel index n associated with the green line. The same applies to the upper right (valley), the bottom left (wall), and the bottom right (cliff) profiles. 
However, it is not straightforward to use Eq. (4) to control the error ℰ via changing the kernel index n in ℐ(n), owing to the dependence of the expression on the second derivative of the function. Instead, we introduced an estimator parameter λ, such that λ = 1 when the system behaves linearly, but λ> 1 in regions where the fluid departs from linearity. The local value of λ was then used to set the exponent n of the kernel that reduces the error in the density estimation. A similar strategy was introduced by Sigalotti et al. (2006) to select the value of the smoothing parameter h. We comment on the similarities and differences between our proposal and that of Sigalotti and coworkers in the concluding section.
For each particle a an estimator λ_{a} is defined, (8)where , where ρ_{b} is a density estimate calculated with a kernel index n, and is the number of neighbors of the particle. When the values of λ_{a} are known, a new kernel index n_{a} is assigned to each mass point according to (9)where n_{0} is a constant baseline value of the kernel index set at the beginning of the calculation, Δn is the highest allowed jump above the baseline value, and λ_{c} is a scaling parameter (λ_{c} ≃ 0.5). At each model a trial value of h and n are picked and the density is computed. These trial values are iteratively refined, with the scheme explained in Appendix A, until the constraints on h (h^{d}ρ = const.) and n (Eq. (9)) are fulfilled. n = n_{0} can usually be taken everywhere for the first model, but it departs from n_{0} in fluid regions with a nonlinear behavior. Below we refer to the adaptive sinc kernel indexes as n(x) in 1D numerical experiments with static configurations, n(x,t) for timedependent 1D hydrodynamic simulations, and n(r,t) in more than one dimension.
The function f(ξ) must fulfill at least two limiting conditions: (10)This ensures that the kernel index remains close to its baseline value, n_{0}, in regions where but becomes n ≃ n_{0} + Δn in regions with a clearly nonlinear behavior.
A suitable function f(ξ) used in this work is (11)This function has the interesting property that df(ξ) / dξ) ≃ 0 at low ξ, making it insensitive to numerical noise. Nevertheless, the function becomes steep at moderate ξ while flattering again at ξ ≫ 1.
An important feature of this scheme is that it is compatible with the Lagrangian derivation of the SPH equations. In other words, the gradient of the kernel index can be incorporated into the Euler equations in the same way as the gradient of the smoothing length is taken into account in the standard SPH (Springel & Hernquist 2002). The reason is that the estimator λ_{a} defined by Eq. (8) admits an explicit derivative with respect to ρ_{a}, (12)Using Eqs. (9), (11), and (12), is straightforward to compute , which is needed to correct the Euler equations from the new gradn terms (see Sect. 4). The constraint in n set by Eq. (9) is, however, of a different kind than that arising from h^{d}ρ = const. (used to set the value of the smoothing length at each time step). While the latter is a real physical constraint and a direct consequence of mass conservation, the former arises from a mathematic consideration of linearity in a local fluid region.
3.1. Fitting sharp 1D density profiles
In this section we consider the ability of the proposed scheme to reproduce the 1D density profile of some mass distributions that often appear in hydrodynamics calculations. These profiles are referred to as mountain, valley, wall, and cliff (see Fig. 4). Their mathematical expressions are
where δ is the characteristic width of the function. The convolution of this curve with the kernel provides the SPH density values. The density is calculated with the standard SPH summation (17)The parameter values in expressions (13–16) are shown in Table 2, and δ<h, . The results of the calculations are depicted in Fig. 4. These profiles are idealized mathematical curves mimicking physical situations of considerable interest in gas dynamics. A mountainlike profile with density contrast of four can appear in regions with a strong shock moving through a perfect gas with adiabatic index γ = 5 / 3. The inverted Gaussian (valleylike) structures may appear during the propagation rarefaction waves. Walllike structures can be found at fluid regions that contact rigid boundaries. Selfgravitating bodies usually end in rarefied atmospheres with steep (clifflike) density gradients. In all the curves shown in Fig. 4, the characteristic width in the steepest regions is lower than the smoothing length. Thus we expect problems in the SPH approach for particles located in the neighborhood of the discontinuities.
The analysis of these idealized curves unambiguously indicates an enhancement of the numerical fitting when adaptive kernel indexes are used. The improvement is especially good in the lowdensity regions that host steep gradients (see, for example, the upper right panel in Fig. 4 that shows an inverted Gaussian). These rarefacted regions are precisely the regions where the standard SPH gives the poorest results because the smoothing length h becomes longer to satisfy the constraint ρh^{d} = const.
The comparison between the mountain and valleylike profiles suggests that the effect of equalizing the error is not symmetric. In the first case the highest value of n (≃8) is achieved not at the peak of the Gaussian, but at the base of the profile where the curve becomes flat, while in the second case n ≃ 8 is taken just at the bottom of the inverted Gaussian. In both cases the index n clearly increases in regions where the second derivative does not vanish. The walllike case is similar to that of the mountain, but with a plateau to the right of the profile. Again it can be seen as n peaks at the base of the wall where the second derivative is larger, maximizing the differences among the profiles calculated with n = 3, n = 5 and n(x). Finally, the case of the cliff, depicted in the bottomright panel of Fig. 4, is particular because the particle sample near the base of the cliff is usually sparse, if not void (especially in 3D applications). Even though the true profile is not reproduced by any of the calculations, we see that the better fit is achieved using variable exponents. This suggests that using adaptive sinc kernels might be of great interest to simulate phenomena in the envelope of selfgravitating bodies, as long as a sufficient sample of particles is available.
Fig. 5 Profiles of density, velocity, zoom of density and kernel index of the 1D blastwave at time t = 0.08 s. The continuum red line is the analytical value. Pink (dots), blue (dashed), and green (dashed) denote n = 3, n = 5, and n adaptive, respectively. The profile in black (dashed) lines and crosses (density zoom) plots n adaptive, but keeping n = 5 in the artificial viscosity terms. 
4. Hydrodynamic equations
The numerical scheme described above was validated through several hydrodynamic tests and compared with the results obtained from keeping the kernel exponent unchanged. We used the standard SPH written in the Lagrangian formulation as described in Monaghan (2005), Rosswog (2009), Springel (2010b), and Price (2012).
4.1. Euler equations with gradh and gradn terms
Small changes of the standard SPH scheme are necessary to incorporate the sinc family of kernels with adaptive indexes. The Euler equations do not change: (18)(19)(20)where is given by Eq. (2), and the remaining symbols have their standard meaning (Monaghan 2005). The parameter Ω_{a} includes the relevant information to compute not only the gradh, but also the gradn derivatives: (21)where . The last term on the RHS in Eq. (21) accounts for the correction for the gradient of the exponent of the sinc kernel n, which the distinctive feature of our proposal. The derivative (22)can be computed from Eqs. (6), (8), and (9). Details of the implementation of the gradh and gradn corrections are given in Appendix A.
For the artificial viscosity we used the recipe described in Monaghan (1997), inspired by the Riemann solvers formulation, where the term Π_{ab} accounting for the viscous pressure is (23)where is an estimate of the signal velocity between particles a and b and w_{ab} = r_{ab}·v_{ab}/  r_{ab}  is the relative velocity projected onto the separation vector. Following Springel (2010b), we used a constant α = 4 / 3 to carry out the simulations described below, so that Π_{ab} remains close to the classical SPH artificial viscosity introduced by Monaghan & Gingold (1983). This particular form of the artificial viscosity has the advantage that there is no explicit dependence of viscosity on the smoothing length, because using n(r,t) makes h less reliable as an indicator of resolution. In principle, the viscous terms in Eqs. (19), (20) could be computed using a different kernel index than those depending on gas pressure. We have not found relevant differences among the results of the tests described below when the actual index n_{a} of the particle, calculated with expression (9), or the constant baseline value n_{0} is used to compute the viscous part of Euler equations. The only exception was the 1D blast wave test, where the variable exponents n_{a}(x,t) led to a narrower spike in the density peak. For that reason variable kernel indexes were also used to estimate the contribution of viscous terms to momentum and energy.
The calculation of the Euler equations is preceded by a brief preconditioning stage, where the optimal values of h and n are set. The value of h is chosen so that the mass within a volume h^{d} is constant during the calculation. An initial pilot value of the density ρ_{a}, calculated with a trial n_{a}, as well as are evaluated at this point. Then the selfconsistent new values of h_{a} and n_{a} are found using the NewtonRaphson (NR) iterative scheme described in Appendix A. Regardless of setting Δn = 0 in Eq. (9) or imposing , the preconditioning algorithm is restored to the standard description, in which the kernel index is kept constant, and the smoothing length and density are jointly updated.
The value of the freeparameter λ_{c} sets the sensitivity of the kernel index with respect to λ. For [ n_{0},λ_{c},Δn ] in Eq. (9) we used [ 5,0.5,5 ] in all the tests below, which yielded satisfactory results. A lower value of λ_{c} leads to larger exponents, but also increases the noise level.
5. Hydrodynamic tests
5.1. Onedimensional tests
5.1.1. Blast waves
Reproducing a strong 1D blast wave with a known analytical solution is a powerful test for any hydrocode. The main goal here is to analyze if the adaptive kernel index algorithm is robust and leads to results better than or at least comparable with the calculation with constant exponents. Our first test was carried out with the same initial setting as in Monaghan (1997). From now on, initial models are specified by [ ρ,vel,γ,P,Δ ], where γ is the constant that relates pressure and internal energy, P = (γ − 1)ρu, and Δ the interparticle distance. For this test [ 1,0,1.4,10^{3},0.005 ] for x< 0 and [ 1,0,1.4,0.01,0.005 ] for x> 0. Simulations were carried out using constant kernel indexes n = 3 and n = 5, as well as variable kernel indexes. A model was also run with variable exponents n_{a}, but keeping n_{a} = n_{0} in the viscous terms of momentum and energy equations.
The results of the calculations are summarized in Fig. 5, where we show the profiles of density, velocity, and kernel index at t = 0.08 s. There are no substantial differences between the different models. They all agree well with the analytical profile. From the fine details, however, we see that the calculations with n(x,t) depict the density in the rarefaction tail of the wave slightly better (bottom left panel in Fig. 5). In this case, we see a small spike that only affects one particle in the plateau at highest density. This feature disappears if constant n_{0} is taken to compute ∇W_{ab} in the viscous terms of Eqs. (19), (20).
Fig. 6 Density, kernel index, zoom of density and pressure profiles of the 1D shock wave born from a single particle. The details are the same as in Fig. 5. 
The algorithm to selfadapt n(x,t) is robust and works very well, detecting strong gradients of density and interphases, as suggested in the bottom right panel of Fig. 5. The index of the sinc kernels changes only in a very narrow region at the sides of the density peak, almost reaching its highest allowed value n = 10. Note also the similarities with the mountainlike static profile of Fig. 4, where n(x,t) peaks twice around the maximum in the density profile. The density profiles of models with different n are also similar at the lowdensity tail in the shocked region.
As a variation of the previous test, we tracked the evolution of a 1D pointlike explosion. In this case, the density contrast between the peak and the bottom of the profile is higher than in the preceding case. We started from a homogeneous distribution of particles with [ 1,0,1.4,0.01,2 × 10^{3} ]. The explosion was initiated by increasing the internal energy of the central particle by a factor 10^{6}. As before, the evolution was followed using three prescriptions for the kernel index, n = 3,n = 5, and n adaptive. Figure 6 depicts the density profile at time t = 0.016 s for the different indexes and initial resolution h_{0} = 1.5Δ. The pattern consists of two strong shockwaves moving in opposite directions, separated by a diluted region. Again we see that n changes abruptly around the discontinuities. Nevertheless, the density profile matches well, regardless of the value chosen for the kernel index. The exception is the central diluted zone, which is better described when equalization is turned on, as shown in the bottom left panel of Fig. 6. The profile of pressure (normalized to the pressure peak at the shock front) is also depicted in the same figure and compared with the analytical profile. Around the peak of the blast all cases agree well, but this is different at the central, lowpressure region. Still, the calculation with adaptive index provides a better approach to the pressure in that zone. This feature is also seen in two dimensions, as commented in Sect. 5.2.1.
5.1.2. Shocktube test
This is a similar test as before, but now the shock and the rarefaction waves are much weaker. In this case, a box is filled with a gas so that the pressure in the leftmost part of the box is higher than in the right side. At t = 0 s both regions are separated by a wall. When the wall is removed, the two regions begin to mix and a shock wave appears that moves through the lowpressure region, while a rarefaction wave digs into the highpressure zone. The initial conditions are left [ 1,0,1.4,1,2.5 × 10^{4} ], right [ 0.125,0,1.4,0.1,2 × 10^{3} ].
Fig. 7 Characteristic profiles of density, specific internal energy, pressure, and velocity of the 1D shocktube problem during the selfsimilar evolution. The continuum red line is the exact solution, the green line was calculated with n adaptive. 
A summary of the results is given in Fig. 7 where the profiles of density, internal energy, pressure and velocity are shown and compared with the analytical values. In this test the resulting profiles calculated with constant n = 3, n = 5 and n adaptive are nearly identical and only the result for n(x,t) is given. This matches the exact profile very well. Nevertheless, we also see a sharp spike in internal energy and pressure at the contact discontinuity. This feature (also present in the calculations with n = 3, n = 5) is known to show up when the contact discontinuity is not smoothed at t = 0 s and there is no heat diffusion term, driven by an artificial conductivity, included in the energy equation.
In Fig. 8, we show the evolution of the profile of the kernel index n(x,t). The highest values of n are achieved at t ≃ 0 when the density contrast around the contact discontinuity is highest and its profile steep. Nevertheless, they decay fast to values close to the baseline value n_{0} = 5 as soon the selfsimilar state is achieved, which makes the results very similar to those obtained with constant n.
5.1.3. Sjögreen test
As described by Einfeldt, et al. (1991), this gasdynamics problem involves the propagation of two symmetric rarefaction waves through a perfect gas with γ = 1.4. The Sjögreen test can be easily handled with SPH, but not with methods using iterative Riemann solvers unless special techniques are used. To initiate these waves, the initial conditions were set as in Monaghan (1997), with one half of the system moving to the right with [ 1,2,1.4,0.4,0.001 ], while the other half moves to the left with [ 1, − 2,1.4,0.4,0.001 ]. As a result, a cavity filled with a very diluted gas appears at the center of the system. The geometry of this fluid cavity resembles the valleylike profile depicted in the upper left panel of Fig. 4. A comparison between the profiles of several magnitudes at t = 0.9 s, obtained with and without equalization, is provided in Figs. 9 and 10 for two values of the initial smoothing length h_{0} (note that a logarithmic scale was used to highlight the differences in the diluted region). Even though the results are good in all cases, there is a clear improvement when the equalization algorithm is included, especially for h_{0} = 3Δ. The lowest values achieved by the density, internal energy, and pressure are closer to the analytical expectations. These results agree qualitatively with the static valleylike case of Fig. 4. The velocity profile at the center is slightly flatter when the adaptive kernel index is used, being also closer to the analytical solution.
The profiles of h and n for the case h_{0} = 1.5Δ are shown in Fig. 10. Both the smoothing length and the adaptive kernel index steeply increase in the vicinity of the density minimum. The variable index n is thus controlling the loss of resolution caused by the growth of h. This effect is not linear, however, because the model calculated with equalization has a higher value of h at the lowest density than models calculated with constant n. This is a consequence of the strong coupling between ρ,h, and n which, as mentioned above, are selfconsistently found at each step using an iterative NR scheme.
5.2. Multidimensional tests
5.2.1. Sedov test in 2D
To study the evolution of a spherical SedovTaylor blast wave, we conducted a test that involve the propagation of a deltafunction signal. This gasdynamical problem includes a pointlike explosion inside a homogeneous system. The explosion rapidly evolves towards a selfsimilar wave with a known analytical solution (Sedov 1959). This is a very demanding test in more dimensions than 1D, where the resolution is usually too low to yield reliable values of the magnitudes around the peak of the wave or close to the origin of the explosion. We wish to know if the combined adaptive h − n scheme can describe this phenomenon better. To trigger the explosion, a δlike function was imposed on the internal energy at t = 0 s, (24)where r is the distance to the explosion center and u_{0} = 10^{7} erg g^{1}, σ = 0.02 cm. For this test the initial interparticle separation was Δ = 4 × 10^{3}. The initial value of the smoothing length was set to encompass neighbors. The profiles of several magnitudes during the selfsimilar evolution of the blast are shown in Fig. 12. As in the preceding tests, the equalization mostly affects the shocked region, although its imprint is not strong. The rarefacted tail of the blast wave is better described when variable n(r,t) are used. In particular, the pressure profile downstream shows a clear dependence on the index of the sinc kernels. The conservation of energy is quantified in the bottom right panel of Fig. 12 with the adaptive n(r,t) scheme providing the best results. It is interesting to note that although the SPH formalism is built to exactly conserve energy, the conservation is usually not perfect in practice owing to the small errors in particle localization, especially when neighboring particles have very different smoothing lengths. Using a large n in S_{n} makes interpolations less susceptible against small fluctuations at the outer edge of the kernel. In this sense, including the equalization lowers these errors and improves the total energy conservation. The profile of n(r,t), when equalization was included, is depicted in the bottom left panel of Fig. 12 and in Fig. 13. We see two regions where the kernel index becomes higher than its baseline value n_{0} = 5, one around the shock front peaking at n ≃ 6.5 and other at the postshock diluted region with a highest value of n ≃ 9.5.
Fig. 8 Profiles of n at different times for the shocktube numerical experiment. 
Fig. 9 Density, internal energy, pressure, and velocity profiles for the Sjögreen test with initial particle separation h_{0} = 1.5Δ. The details are the same as in Fig. 5. 
Fig. 11 Profile of the smoothinglength, h(x), for n = 3 (dashed pink line), n = 5 (long dashed blue line) and n adaptive (continuum green line) for the Sjögreen test. The profile of n(x) is also shown (black dots). 
Finally, a calculation was launched with a high value of the kernel exponent, n = 10, in all particles. A color density map for cases n = 3, n = 5, n = 10 and n adaptive is provided in Fig. 14. For n = 10 the density distribution is affected by the initial particle setting in a rectangular lattice (sometimes referred to as hourglass instability). As expected, the spherical symmetry is better preserved for the loworder interpolator n = 3, but the cases n = 5 and n adaptive are also very good. We conclude that the use of high kernel indexes must be reserved to fluid regions that host sharp density gradients (see also Sect. 5.2.2). Loworder interpolators are more efficient in suppressing numerical noise, but they are less accurate and more prone to undergo pairing instability. A conservative procedure is to take a moderate index n in all fluid regions to reduce the numerical noise, but switch to a larger n wherever a discontinuity is found.
5.2.2. KelvinHelmholtz instability in 2D
The KelvinHelmholtz instability appears when there is a sufficient shear velocity in the interface layer between two fluids with different densities. Small perturbations of the velocity field in the orthogonal direction to the interface emerge and lead to a mixing of the two fluids. This is usually simulated in a box with periodic boundary conditions, where two fluid regions are defined with densities ρ_{1} and ρ_{2}. The two layers have opposite parallel velocities, which leads to a shear discontinuity in the contact interface. To develop the instability, a small perturbation is seeded at the interface as a sinusoidal mode of length scale L. Recent SPH simulations of the KH instability can be found in McNally et al. (2012) and Hopkins (2012).
We simulated a central band of a highdensity fluid ρ_{1} moving in a lowdensity medium ρ_{2} in a squared lattice of 1 cm side in the XY plane using N = 62,500. The density around the interface was not smoothed. The initial setting was [ 1, − 0.5,5 / 3,2.5,0.006 ] for y ≤ 0.25,y ≥ 0.75 and [ 2, + 0.5,5 / 3,2.5,0.003 ] for 0.25 <y< 0.75. The initial smoothinglength was chosen so that every particle sees neighbors.
A sinusoidal perturbation of the v_{y} component of the velocity field was seeded at t = 0. Then, for the initial velocity we have (25)where we took m = 2 and Δv_{y} = 0.01 cm s^{1}, a small perturbation indeed.
First of all, we would like to stress that the calculation with n = 3 was a complete failure because the perturbation failed to emerge. The reason for this was that according to Fig. 3, the initial number of neighbors is much higher than necessary to suppress pairing instability. In the calculations with n = 5, however, the pair , lies only moderately above that line, and no trace of particle clustering was detected during the simulation. Particle clustering can also be avoided, even for n ≃ 3, using a a different SPH approach to the fluid equations, such as those based on an integral approach to the derivatives (IAD; GarcíaSenz et al. 2012).
Fig. 12 Profiles of ρ,P,n (averaged in concentric shells) and evolution of energy conservation during the 2D blast wave propagation (Sedov test). The continuum line in red is the classical Sedov solution, and profiles in pink (dots), blue (dashed) and green (longdashed) are for cases n = 3, n = 5 and nadaptive, respectively. 
Figure 15 shows a density color map of the growth of the KelvinHelmholtz instability at different times for the calculations using n = 5, n = 10 and n adaptive. The simulation with n = 10 is manifestly poorer, suggesting again (see the preceding section) that choosing a large n from the beginning is not a safe option. Cases n = 5 and n are almost indistinguishable; both lead to a clear growth of the instability with its characteristic pattern. The last row of Fig. 15 shows a color map of n(r,t) for the same models as depicted in row 3. The tracking algorithm perform well because n only increases in a thin shell around the interface. Nevertheless, the change in n is not large and the dynamical evolution remained close to that with n = 5.
5.2.3. Astrophysical application: gravitational collapse of a polytrope
Finally, we simulated the gravitational collapse of a Sunlike polytrope with and without the equalization algorithm and compared the results with the output of a wellknown 1D Lagrangian hydrocode (the AGILE hydrocode by Liebendörfer et al. 2002). We carried out the 1D models with AGILE taking 260 grid points. This provides an output with much better resolution than that of the SPH hydrocode and serve as a suitable reference model.
The initial model used for the comparison was a 1 M_{⊙} spherically symmetric polytrope of index 3. The radius was set to 1 R_{⊙}, which results in a central density of ρ_{c} = 76 g cm^{3}. We built specific equilibrium initial models for each case (with and without equalization) by distributing N = 10^{5} particles in 3D according to the 1D density profile, and let them relax to the hydrostatic equilibrium. The EOS used in the simulations was that of a perfect gas with γ = 5 / 3. The initial value of h in the SPH calculations was chosen to encompass neighbors.
The equilibrium structure was then suddenly destabilized by removing 20% of its internal energy, so that the star collapsed under the force of gravity. At some point, the collapse in the central zone was halted because of the increase of pressure, and an accretion shock formed that moved through the infalling material to ultimately eject the surface layers of the polytrope. That scenario contains several pieces of physics of great interest because accretion shocks and pulsational instabilities are very common in astrophysics.
The evolution of the central density during the collapse and the rebound of the star is shown in Fig. 16, the profile of several variables at three elapsed times is provided in Fig. 17. In general, all calculated models show a similar behavior during the implosion and first oscillation of the polytrope. The first peak of central density is achieved after ≃15 min in all simulations. Nevertheless, the exact value of the peak is affected because the resolution is higher in the AGILE 1D calculation and lower for SPH with n = 3, as expected. The calculations with n = 5 and n adaptive virtually led to the same maximum in the central density. The discrepancy between n adaptive and the reference 1D model is ≃ 6%. Figure 16 also depicts the evolution of the fraction of total energy lost during the first hour, which remains below 0.2% for all the SPH models. As in the Sedov test, the calculations with higher exponents, n = 5,n(r,t) conserve the energy better than that with the cubicsplinelike kernel, case n = 3.
The profiles of density, velocity, and internal energy at times t = 870 s, 1086 s, and 1311 s are depicted in Fig. 17. The density profiles do not show any significant difference between the n = 5 and n(r,t) calculations. In both cases the discontinuity at the accretion shock is smoothed in a similar way and is less pronounced than in the reference model. The radial velocity profiles are shown in the upper right panel of Fig. 17. They show some differences at the position of the accretion shock; the simulation using n(r,t) better matches the AGILE results at t = 1086 s and t = 1311 s. In particular, the lowest velocity is much better captured when the equalization is included. A similar behavior is observed in the profile of the specific internal energy depicted in the bottom left panel of Fig. 17, but there the differences are not as accentuated as in the velocity profile. The bottom right panel of the same figure shows the distribution of n along the star. The algorithm detects both the accretion shock and the surface of the polytrope at t = 870 s, while for longer elapsed times n(r,t) follows a walllike profile with the baseline value n_{0} = 5 until ≃0.4 R_{⊙} and n ≃ 10 at the surface.
Fig. 13 Rendering of the sinc kernel index n(x,y) for the selfsimilar wave shown in the bottom left panel of Fig. 12. 
Fig. 14 Density color map of the Sedov wave for n = 3 (upper left), n = 5 (upper right), n = 10 (bottom left), and nadaptive (bottom right). The spherical symmetry is poorly preserved for a large constant exponent such as n = 10 in the sinc kernel. 
6. Discussion and conclusions
In the standard formulation of the SPH method the resolution is bounded to the local density value, meaning that the rarefacted zones of the fluid are intrinsically handled with a lower resolution than the highdensity regions. We proposed a method to equalize the error in diluted regions that is robust and easy to implement, with a low computational overload. The formulation of the method relies on the definition of a local estimator of the linearity of density. According to Eq. (8), the definition of that estimator, λ, is fairly simple and its value is used to control the accuracy of the interpolations. A similar method was proposed by Sigalotti et al. (2006) as a way to set the value of the smoothing length at each steptime. Our proposal differs from that of these authors in several ways. First, in our method h is set in the standard manner, keeping the mass constant around a particle, while the value of λ sets the value of the exponent, n, of the sinc kernels. Second, unlike Sigalotti et al. (2006), who neglected the gradh terms, we included the gradh and gradn corrections to the momentum and energy equations. The computation of these corrections is compatible with the Lagrangian derivation of the fluid equations. Third, the specific mathematical expressions used to set n are different from those used by Sigalotti and coworkers to set the value of h. In our proposal, we constrained n to the range n_{0} ≤ n ≤ n_{0} + Δn, with the boundaries achieved asymptotically.
Fig. 15 Density color map showing the growth of the KelvinHelmholtz instability for cases n = 5 (first row), n = 10 (second row), n adaptive (third row). The last row plots the contours of n(x,y) for the adaptive case from n = 5 in the black zones to around n = 7 in the brightest red zones. 
The proposed algorithm works well with static 1D particle distributions. According to Fig. 4, the zones with sharp density gradients are better described using an adaptive kernel index n(x). Unlike the adaptive h(x), the improvement due to n(x) is more pronounced in the lowdensity tail of the profiles. Therefore using both h(x) and n(x) tends to equalize the error along the system.
In hydrodynamic calculations, however, the improvements are not as pronounced as in the static profiles. The main reason is that sharp density gradients are smoothed by the artificial viscosity that widens the discontinuities to twice or thrice the smoothing length. Moreover, the mechanism by which h and n selfadapt is not longer linear, and sometimes a selfconsistent increase of n is followed by an increase in the number of neighbors, making the enhancement in resolution less noticeable (see, for instance, Fig. 11). Still, the hydrodynamic tests confirm the main results attained with 1D static profiles: a moderate improvement in the description of the rarefacted regions of the gas, usually attached at the rear tail of shock waves. The careful handling of these postshock regions must not be disregarded because it is as important as the shock front itself: in these tails hydrodynamic instabilities may grow under the appropriate physical conditions (for example, the RayleighTaylor instability in the regions between the forward and reverse shocks in supernova remnants).
The search for the optimal n(r,t) can be made in the same NR loop as was used to update h(r,t) with very little changes. At each iteration the value of and the local arithmetic mean of lnρ have to be stored, but the overload is small if a list of the neighbors of each particle is stored in an array and used to localize particles when necessary. A switch can be used to include or exclude the equalizing option, as shown in Fig. A.1. The algorithm is very efficient in detecting discontinuities. It was able to track the contours of shocks, walls, and surfaces in all the tests. The equalization does not interfere with the development of the KelvinHelmholtz instability either because it neither enhances nor diminishes the growth.
Fig. 16 Trajectory of the central density during the implosion of a Sunlike polytrope and percent of energy conservation during the first elapsed hour. The red continuum line is the 1D calculation with the implicit Lagrangian hydrocode AGILE. Dashed lines are for n = 3, n = 5, and n adaptive. Light continuum lines in pink, blue and green show the percent of energy conservation. 
Fig. 17 Profiles of density, velocity (in units of 10^{7} cm s^{1}), internal energy, and index of the sinc kernel at times t_{1} = 870 s, t_{2} = 1086 s, t_{3} = 1311 s, corresponding to the collapse of a Sunlike polytrope. The lines in green, pink, and orange are for n adaptive and the blue, light blue, and black lines for n = 5 at the same elapsed times. The continuum lines in red have been obtained using a Lagrangian hydrocode with spherical symmetry. 
The application to a specific astrophysical problem: the collapse and subsequent rebound of a Sunlike polytrope was also satisfactory. The calculation with equalization led to better profiles of velocity and specific internal energy with an adaptive n(r,t) increasing in the rear of the accretionshock front and at the surface. Nevertheless, while a high value of n at the shock is driving a clear enhancement of the internal energy and radial velocity profiles, its impact on the surface layers was weak. The reason is that interpolations at the boundaries of selfgravitating bodies are not as accurate as in the interior because of the scarcity of sampling points in the outermost regions of the envelope. To adequately solve the surface layers in 3D and estimate the real effect of equalization a huge increase in the number of particles would be necessary.
Among other advantages, the sinc family of interpolators introduces an additional degree of freedom to control the resolution in SPH. The simultaneous (implicit or semiimplicit) search for h and n increases the computational burden, but this is no great concern unless very many particles require a hard refining of n. In this respect, we estimated a ≃10% overload in the simulation of the 2D Sedov pointlike explosion. The computational penalty will be weaker in current astrophysical scenarios where gravity and/or a complex physics are incorporated in the numerical scheme.
A priori, working with n(r,t) can also be a potential source of numerical noise, which may affect the development of small fluid instabilities. In this respect, we found no spurious effect in the growth of the KH instability, but more work is needed to confirm this last point. Additionally, other functional forms of the estimator λ_{a}, different from that used in this work given by Eq. (8), might be devised to control n and better adapt the abilities of current SPH codes to handle specific physical problems.
Acknowledgments
This work has been funded by the Spanish MEC grants AYA201015685, AYA201123102 and DURSI of the Generalitat of Catalunya (D.G.S. and J.A.E.). R.M.C. acknowledges the support by the Swiss Platform for HighPerformance and HighProductivity Computing (HP2C) within the supernova project and the Platform for Advanced Scientific Computation (PASC) within the DIAPHANE project. R.M.C. and K.E. were also supported by the ERC grant FISH. D.G.S. was also supported by the EuroGENESIS and CompStar progams. The rendered SPH plots were made using the freely available SPLASH code (Price 2007).
References
 Cabezón, R. M., GarcíaSenz, D., & Relaño, A. 2008, J. Comput. Phys., 227, 8523 [NASA ADS] [CrossRef] [Google Scholar]
 Cabezón, R. M., GarcíaSenz, D., & Escartín, J. A. 2012, A&A, 545, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068 [NASA ADS] [CrossRef] [Google Scholar]
 Einfeldt, B., Munz, C. D., Roe, P. L., & Sjögreen, B. 1991, J. Comp. Phys., 92, 273 [Google Scholar]
 GarcíaSenz, D., Cabezón, R. M., & Escartín, J. A. 2012, A&A, 538, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375 [Google Scholar]
 Hopkins, P. E. 2013, MNRAS, 428, 2840 [NASA ADS] [CrossRef] [Google Scholar]
 Liebendörfer, M., Rosswog, S., & Thielemann, F. K. 2002, ApJS, 141, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 McNally, C., Lyra, W., & Passy, J.C. 2012, ApJ, 201, 18 [NASA ADS] [Google Scholar]
 Monaghan, J. J. 1992, ARA&A, 365, 199 [Google Scholar]
 Monaghan, J. J. 1997, J. Comput. Phys., 136, 298 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Monaghan, J. J. 2005, Rep. Prog. Phys., 68, 1703 [Google Scholar]
 Monaghan, J. J., & Gingold, R. A. 1983, J. Comput. Phys., 52, 374 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. 2007, PASA, 24, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. 2012, J. Comput. Phys., 231, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Rosswog, S. 2009, New Astron. Rev., 53, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Rosswog, S. 2014 [arXiv:1405.6034] [Google Scholar]
 Saitoh, T., & Makino, J. 2013, ApJ, 768, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academic Press) [Google Scholar]
 Sigalotti, L. D. G., López, H., Donoso, A., Sira, E., & Klapp, J. 2006, J. Comput. Phys., 212, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2010a, MNRAS, 401, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2010b, ARA&A, 48, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649 [NASA ADS] [CrossRef] [Google Scholar]
 Valdarnini, R. 2012, A&A, 546, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wendland, H. 1995, Adv. Comput. Math., 4, 389 [CrossRef] [MathSciNet] [Google Scholar]
Appendix A: Implementation of the algorithm to compute gradh and gradn
A suitable mathematical expression giving the corrections by the gradh and gradn terms can be obtained as a simple extension of the reasoning used to compute the gradh terms (Rosswog 2009). The discretized fluid movement equations are derived using the EulerLagrange formulation (A.1)where r_{a} and v_{a} refer to the position and velocity of particle a. The Lagrange function of the system is (A.2)where u_{b},s_{b} are the specific internal energy and entropy of particle b. Inserting Eqs. (A.2) into (A.1) and admitting isentropic evolution, the movement equations for particle a are written (A.3)Following Rosswog (2009), the density gradient in Eq. (A.3) (also needed to compute the energy equation) is calculated as (A.4)where ∇_{a} is the derivative with respect to the spatial coordinates and (A.5)with and (A.6)which can be computed from Eqs. (9), (11) and (12). Inserting Eqs. (A.5) and (A.4) into (A.3), the form of the momentum equation used in this work, Eq. (19), is easily recovered.
Fig. A.1 Flow chart of the preconditioning algorithm built to set the values of ρ,h,n before starting the hydrodynamics. 
The ensuing algorithm to update h,n and calculate the gradh, gradn corrections of particle b was implemented using an NR iterative scheme, (A.7)
where is a constant, set at the beginning of the simulation, and f(ξ_{b}) is defined in Eq. (11), hence, (A.8)Note that for Δn = 0 the NR reduces to the standard scheme, where the kernel exponent is kept constant and h is updated according to the local density. For [ n_{0},Δn,λ_{c} ] we have taken [ 5,5,0.5 ], which led to reliable results in the numerical tests with neighbors in 2D and 3D, respectively. To speed up the calculations, it is highly recommended to store the values of (the ln [ sinc ] can be used to compute ) as a function 0 ≤ v ≤ 2 in an array, and interpolate from them to obtain any kernelrelated magnitude. A sample of 2 × 10^{4} equally spaced points was good enough for all tests presented in this work.
A flow chart of the preconditioning moduli is given in Fig. A.1. When implementing the algorithm, particles that have already converged need to be carefully removed from the general NR loop, taking them into account only to compute . If the algorithm is well balanced and optimized, the computational overload should remain at a few %, unless very many particles require hard refining.
All Tables
All Figures
Fig. 1 Left: profile of several sinc (S_{n}, continuum lines) and Wendland kernels (ψ_{l:k}) (rescaled to a common range [0, 1]). Points + (in red), × (in green) and ∗ (in blue) are for ψ_{3:1}, ψ_{4:2} and ψ_{5:3} respectively. Right: same as before, but for sinc and cubic, quartic, and quintic spline kernels. 

In the text 
Fig. 2 Fourier transform,  ℱ_{3}  of sinc kernels with n = 3,4 and 10 as well as the cubic spline, where the dashed lines indicate the negative portions of the curves. They can be compared with the Fourier transform of several Wendland and spline kernels given in Fig. 2 by Dehnen & Aly (2012). 

In the text 
Fig. 3 Color map of the logarithm of the error in density estimation in 2D (left) and 3D (right) as a function of the kernel index and number of neighbors (see Sect. 3 for a complete explanation of this diagram). The dashed line in red is the rough critical limit separating the region susceptible to particle pairing (above the line). The blue dotted line denotes the region where the approach of integrals by summations in density calculation is too sensitive to particle distribution (below the line). 

In the text 
Fig. 4 Fitting of several 1D density profiles using the sinc family of kernels with constant and selfadaptive indexes. The upper left panel shows a Gaussian, mountainlike, profile. The exact (e) analytical value is shown in red, the result with the variable (v) kernel index, constant n = 5 and n = 3 in green, blue, and pink, respectively. The light blue line shows the profile of the smoothinglength normalized to its highest value (achieved just at the limits of the system, x = 0 and x = 1). The orange line shows the profile of the kernel index n associated with the green line. The same applies to the upper right (valley), the bottom left (wall), and the bottom right (cliff) profiles. 

In the text 
Fig. 5 Profiles of density, velocity, zoom of density and kernel index of the 1D blastwave at time t = 0.08 s. The continuum red line is the analytical value. Pink (dots), blue (dashed), and green (dashed) denote n = 3, n = 5, and n adaptive, respectively. The profile in black (dashed) lines and crosses (density zoom) plots n adaptive, but keeping n = 5 in the artificial viscosity terms. 

In the text 
Fig. 6 Density, kernel index, zoom of density and pressure profiles of the 1D shock wave born from a single particle. The details are the same as in Fig. 5. 

In the text 
Fig. 7 Characteristic profiles of density, specific internal energy, pressure, and velocity of the 1D shocktube problem during the selfsimilar evolution. The continuum red line is the exact solution, the green line was calculated with n adaptive. 

In the text 
Fig. 8 Profiles of n at different times for the shocktube numerical experiment. 

In the text 
Fig. 9 Density, internal energy, pressure, and velocity profiles for the Sjögreen test with initial particle separation h_{0} = 1.5Δ. The details are the same as in Fig. 5. 

In the text 
Fig. 10 Same as in Fig. 9, but for initial particle separation h_{0} = 3Δ. 

In the text 
Fig. 11 Profile of the smoothinglength, h(x), for n = 3 (dashed pink line), n = 5 (long dashed blue line) and n adaptive (continuum green line) for the Sjögreen test. The profile of n(x) is also shown (black dots). 

In the text 
Fig. 12 Profiles of ρ,P,n (averaged in concentric shells) and evolution of energy conservation during the 2D blast wave propagation (Sedov test). The continuum line in red is the classical Sedov solution, and profiles in pink (dots), blue (dashed) and green (longdashed) are for cases n = 3, n = 5 and nadaptive, respectively. 

In the text 
Fig. 13 Rendering of the sinc kernel index n(x,y) for the selfsimilar wave shown in the bottom left panel of Fig. 12. 

In the text 
Fig. 14 Density color map of the Sedov wave for n = 3 (upper left), n = 5 (upper right), n = 10 (bottom left), and nadaptive (bottom right). The spherical symmetry is poorly preserved for a large constant exponent such as n = 10 in the sinc kernel. 

In the text 
Fig. 15 Density color map showing the growth of the KelvinHelmholtz instability for cases n = 5 (first row), n = 10 (second row), n adaptive (third row). The last row plots the contours of n(x,y) for the adaptive case from n = 5 in the black zones to around n = 7 in the brightest red zones. 

In the text 
Fig. 16 Trajectory of the central density during the implosion of a Sunlike polytrope and percent of energy conservation during the first elapsed hour. The red continuum line is the 1D calculation with the implicit Lagrangian hydrocode AGILE. Dashed lines are for n = 3, n = 5, and n adaptive. Light continuum lines in pink, blue and green show the percent of energy conservation. 

In the text 
Fig. 17 Profiles of density, velocity (in units of 10^{7} cm s^{1}), internal energy, and index of the sinc kernel at times t_{1} = 870 s, t_{2} = 1086 s, t_{3} = 1311 s, corresponding to the collapse of a Sunlike polytrope. The lines in green, pink, and orange are for n adaptive and the blue, light blue, and black lines for n = 5 at the same elapsed times. The continuum lines in red have been obtained using a Lagrangian hydrocode with spherical symmetry. 

In the text 
Fig. A.1 Flow chart of the preconditioning algorithm built to set the values of ρ,h,n before starting the hydrodynamics. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.