© ESO 2022

1. Introduction

The smoothed particle hydrodynamics (SPH) is a firmly settled numerical technique that is capable of successfully simulating many cutting-edge problems in physics, astrophysics, and engineering. This technique has undergone a sustained enhancement since its original formulation (Lucy 1977; Gingold & Monaghan 1977) and it is still evolving at a suitable pace (e.g., Springel 2010; Monaghan 2012; Rosswog 2015a, 2020b; Wadsley et al. 2017; Cabezón et al. 2017; Price et al. 2018). A landmark in that evolution concerns the estimation of derivatives and gradients, which can be done using a number of different approaches.

The standard way of calculating gradients is by directly taking the analytic derivative of the interpolating kernel function, which leads to E0-errors, even in the presence of constant pressure fields in non-uniform particle distributions (Agertz et al. 2007; Zhu et al. 2015). One alternative adapts the moving-least squares technique (MLS) to SPH (Dilts 1999) to ensure an exact interpolation of linear functions. In the MLS methods, the SPH kernel is replaced by a new interpolant that is built by combining the original kernel with a linear operator. Such techniques ensure the linear consistency of the interpolations but they also require the solution of a linear system of d + 1 equations (where d is the dimension) for each particle at each integration step. Alternatively, adding renormalization corrections to both the kernel and the gradient of the kernel have been shown to enhance the accuracy in the calculation of gradients and to reduce the tensile instability (Bonet 1999). Tensile instability often arises at boundaries in the flow, separating fluid regions of varying densities. The inconsistency in the treatment of gradients induces an nonphysical surface tension of numerical origin, which inhibits fluid mixing and suppresses hydrodynamic instabilities. A variational principle applied jointly with a kernel renormalization was used by Oger et al. (2007) to reduce the tensile instability and simulate fluid systems with free surfaces. Nevertheless, these MLS and renormalization techniques, in general, do not guarantee the perfect conservation of the whole set of physical laws governing the motion of the fluid, which are at the foundations of the SPH technique. A recent proposal has been the conservative reproducing kernel SPH method (CRKSPH) from Frontiere et al. (2017) and although it is not derived from a Lagrangian, it still enforces a perfect linear interpolation and retains the linear momentum and energy conservation properties.

Another alternative way of estimating gradients was devised by García-Senz et al. (2012). In their proposal, gradients are calculated from an integral expression, so that there is no need to explicitly calculate the analytic derivative of the kernel function. These authors also proved that such an integral approach can be completely compatible with the Lagrangian formulation of the SPH equations, leading to the integral smoothed particle hydrodynamics (ISPH) scheme, named IAD₀ in the seminal paper by García-Senz et al. (2012). It has been shown that the ISPH formulation has the same conservation properties as the standard Lagrangian-derived SPH (Springel & Hernquist 2002). We note that a particularly remarkable feature of ISPH is that it effectively reduces the E0-error in the derivatives (García-Senz et al. 2012; Rosswog 2015a; Valdarnini 2016).

In this work we dig further into the conditions that the ISPH scheme should meet in order to improve the calculation of gradients, which can approach greater precision for linear functions. We found that these conditions are connected with another basic SPH requirement, namely: the correct partition of unit volume. By using one-dimensional numerical experiments, we clarify the link between these two basic properties, namely: any enhancement in the partition of unity leads to a better gradient estimation within the ISPH framework. The results of these 1D toy models are confirmed by detailed 3D hydrodynamic simulations of explosions, collisions, and instabilities.

Here, we propose and discuss a new type of volume elements (VEs) that lead to a better partition of the unity. The ISPH equations of density, movement, and energy are consequently re-formulated, within the Lagrangian framework, so that they become fully compatible with the particular choice of the VEs. The resulting scheme is not only fully conservative, but it also enhances the density estimation and the gradient of linear functions with practically no computational overload.

Thus, the present work is a natural extension of the discussion in Cabezón et al. (2017; hereafter Paper I), where the code SPHYNX was introduced and verified. Here, we focus on modifying the proposal in Paper I, so that the new SPH scheme (Appendix A) incorporates an implicit treatment of the VEs, which is now fully compatible with the Lagrangian derivation of SPH equations, with a consistent treatment of grad-h terms. A second improvement is the use of a novel self-adaptive scheme to surgically reduce the tensile instability in regions with large density gradients. Unfortunately, suppressing that instability requires us to abandon the Lagrangian formulation (Read et al. 2010; Wadsley et al. 2017; Wissing & Shen 2020). In our proposal, the departure from such a formulation is minimal and is limited to regions hosting significant density contrasts. Additionally, we show that our formulation leads to a better behavior in terms of entropy in shocks than other schemes that suppress the tensile instability.

In Sect. 2, we review the main features of the ISPH scheme. We discuss the choice of the generalized volume elements used to compute the summations in Sect. 3. The link between the performance of the ISPH calculation of gradients and the adequate choice of the VEs is highlighted in Sect. 4. Section 5 presents the resulting ISPH equations, while in Sect. 6, we apply our code to several standard tests calculated in three dimensions and we analyze the results in the context of the VEs choice. A summary of our findings and prospects for future works is given in the conclusions section.

2. ISPH formulation

The classical way of evaluating gradients in SPH takes the multidimensional derivative of a function f as:

$$\begin{array}{c}\hfill \mathbf{\nabla }f={\displaystyle {\int }_V}f({\mathit{r}}^{\prime })\mathrm{\nabla }W(|{\mathit{r}}^{\prime }-\mathit{r}|,h)\mathrm{d}{r}^{\prime 3},\end{array}$$(1)

where W(|r′−r|,h) is commonly a Dirac δ-like function named interpolating kernel, which is continuous and derivable, and h is the smoothing length. In the classical SPH formulation, the gradient of the function is estimated by: (1) approaching the integral of Eq. (1) by summations; (2) taking the analytic derivative of the kernel; and (3) assuming that the volume element dr^′3 is adequately represented by m/ρ. Specifically¹,

$$\begin{array}{c}\hfill {⟨\mathbf{\nabla }f⟩}_a={\displaystyle \sum _b\frac{m_b}{{\rho }_b}{f}_b\mathrm{\nabla }W(|{\mathit{r}}_b-{\mathit{r}}_a|,h_a),}\end{array}$$(2)

where a and b refer to neighbouring particles with masses m_a and m_b, respectively. Alternatively, in ISPH, the gradient is calculated from an integral approach (IA), which does not require the explicit analytic derivative of W(|r − r′|, h). A vector integral I(r) is defined as (García-Senz et al. 2012):

$$\begin{array}{c}\hfill I(\mathit{r})\equiv {\displaystyle {\int }_V}[f(\mathit{r}\prime )-f(\mathit{r})]({\mathit{r}}^{\prime }-\mathit{r})W(|{\mathit{r}}^{\prime }-\mathit{r}|,h)\mathrm{d}{r}^{\prime 3}.\end{array}$$(3)

The integral I(r) can be used to find the gradient of a function f(r) in a similar way that the Laplace operator is usually approached from another integral expression in standard SPH (Brookshaw 1985; Monaghan 2005). The IA interpretation of SPH is the consequence of approaching Eq. (3) with summations, along with approaching the function f(r) by a Taylor expansion around the evaluated point,

$$\begin{array}{c}\hfill f({\mathit{r}}_b)-f({\mathit{r}}_a)\simeq \mathrm{\nabla }{\mathit{f}}_a·({\mathit{r}}_b-{\mathit{r}}_a),\end{array}$$(4)

with ∇f_a defined below. The RHS in the integral expression in Eq. (3) becomes:

$$\begin{array}{c}\hfill I({\mathit{r}}_a)={\displaystyle [\sum _b\frac{m_b}{{\rho }_b}f({\mathit{r}}_b)({\mathit{r}}_b-{\mathit{r}}_a)W(|{\mathit{r}}_b-{\mathit{r}}_a|,h_a)]}-\left[f({\mathit{r}}_a){⟨\mathbf{\Delta }\mathit{r}⟩}_a\right],\end{array}$$(5)

where

$$\begin{array}{c}\hfill {⟨\mathbf{\Delta }\mathit{r}⟩}_a={\displaystyle \sum _b\frac{m_b}{{\rho }_b}({\mathit{r}}_b-{\mathit{r}}_a)W(|{\mathit{r}}_b-{\mathit{r}}_a|,h_a).}\end{array}$$(6)

Setting Eqs. (4) and (5) into Eq. (3) allows us to obtain the gradient ∇f_a,

$$\begin{array}{c}\hfill {\left[\begin{array}{c}\partial f/\partial x_1\\ \partial f/\partial x_2\\ \partial f/\partial x_3\end{array}\right]}_a={\left[\begin{array}{ccc}{\tau }_{11}& {\tau }_{12}& {\tau }_{13}\\ {\tau }_{21}& {\tau }_{22}& {\tau }_{23}\\ {\tau }_{31}& {\tau }_{32}& {\tau }_{33}\end{array}\right]}^{-1}\left[\begin{array}{c}{I}_1\\ I_2\\ I_3\end{array}\right],\end{array}$$(7)

where

$$\begin{array}{c}\hfill {\tau }_{ij,a}={\displaystyle \sum _b\frac{m_b}{{\rho }_b}(x_{i,b}-x_{i,a})(x_{j,b}-x_{j,a})W_{\mathit{ab}}(h_a);i,j=1,3,}\end{array}$$(8)

and where the spatial coordinates are represented by x with sub-indexes i or j. Henceforth, we note W_ab(h_a)≡W(|r_b − r_a|,h_a) for the sake of clarity.

It was shown in García-Senz et al. (2012) that Eq. (7) leads to a perfect linear interpolation. Unfortunately, the price to pay is the loss of the full conservation of linear and angular momentum. The exact conservation of the SPH Euler equations and a perfect linear interpolation can only be retrieved simultaneously when ⟨Δr⟩_a → 0. We refer to ISPH as the conservative scheme which simply neglects the term f(r_a)⟨Δr⟩_a in Eq. (5). This is justified because ⟨Δr⟩ is, in fact, the indefinite integral of an odd function, which is zero².

On the other hand, the complete integral approach, which takes into account the f(r_a)⟨Δr⟩_a term in the RHS of Eq. (5), leads to a perfect linear interpolation but is not fully conservative. We refer to it as non-conservative ISPH (ncISPH, hereafter). As commented above, both schemes, ISPH and ncISPH, converge to the same outcome when ⟨Δr⟩_a ≃ 0. Having both, a perfect partition of unity and ⟨Δr⟩_a = 0, has been identified for a long time as the main constraint in order to ensure complete linear consistency in SPH (Liu & Liu 2006).

3. Choice of the volume elements

A common choice with regard to the volume elements in SPH is V_a = m_a/ρ_a, which leads to the widely used density equation,

$$\begin{array}{c}\hfill {\rho }_a={\displaystyle \sum _{b=1}^{n_b}{V}_b{\rho }_b{W}_{\mathit{ab}}(h_a)=\sum _{b=1}^{n_b}{m}_b{W}_{\mathit{ab}}(h_a),}\end{array}$$(9)

which works well provided that the density does not change very much within the kernel range. Nevertheless, the density may appreciably be miscalculated in the presence of shocks and density discontinuities. In these cases, the partition of the unity condition is not fully satisfied:

$$\begin{array}{c}\hfill {\displaystyle \sum _{b=1}^{n_b}\frac{m_b}{{\rho }_b}{W}_{\mathit{ab}}(h_a)\ne 1.}\end{array}$$(10)

The errors in the normalization constraint would introduce a level of uncontrolled errors in the remaining SPH equations. To reduce the normalization errors, the most obvious recipe is to renormalize the kernel itself:

$$\begin{array}{c}\hfill {\rho }_a={\displaystyle \sum _b{m}_b\left(\frac{W_{\mathit{ab}}(h_a)}{{\sum }_c\frac{m_c}{{\rho }_c^0}{W}_{\mathit{ac}}(h_a)}\right),}\end{array}$$(11)

where ${\rho }_c^0={\sum }_b{m}_b{W}_{cb}(h_c)$ is the standard density. A more clever, albeit more complex, variation of that scheme was developed by Colagrossi & Landrini (2003) within the MLS approach, which exactly reproduces the linear variation of a density field (see Gomez-Gesteira et al. 2010 for a review of these topics). Nevertheless, none of those SPH schemes are totally compatible with the Lagrangian formulation of SPH. The Lagrangian formulation of the equations of movement has the advantage that the grad-h terms can be consistently incorporated to the scheme, so that a complete preservation of mass, linear and angular momentum, energy, and entropy is guaranteed. Most MLS applications belong to the realm of computational fluid dynamics (CFD), which usually works with incompressible or weakly-compressible fluids. Having an almost constant density field for the entire simulation implies that the smoothing length also remains constant and the grad-h corrections are negligible (except at the boundaries, which are usually handled with special techniques). For that reason, the MLS methods are mostly used in CFD simulations; however, in the case of astrophysical scenarios, where density contrasts of orders of magnitude are not rare, the Lagrangian approach with a self-consistent treatment of the grad-h terms is preferable. Nevertheless, due to its properties, we are convinced that the methods presented in this paper can be also of use for CFD simulations.

Other options for the volume elements, V_a, may be of interest to address specific problems. The code SPHYNX³ from Cabezón et al. (2017), makes use of the concept of generalized volume elements (Hopkins 2013; Saitoh & Makino 2013). First, a scalar estimator X_a is defined so that the particle volume is:

$$\begin{array}{c}\hfill V_a=\frac{X_a}{k_a},\end{array}$$(12)

with k_a = ∑_bX_bW_ab(h_a). The density of the particle is then calculated as ρ_a = m_a/V_a. Current choices for the estimator that can be found in the literature are $X_a=1,m_a,P_a^k$, where P is the pressure and k ≤ 1. There is, however, a particular choice which, according to Paper I, provides a better partition of the unity, namely:

$$\begin{array}{c}\hfill X_a={\left(\frac{m_a}{{\rho }_a}\right)}^p.\end{array}$$(13)

Setting p = 0 produces the standard volume element for particles with identical mass, whereas 0 < p ≤ 1 gradually improves the kernel normalization when p → 1.

There are two ways to implement the estimator in Eq. (13). The first is to make use of the density calculated in the previous time-step, ${\rho }_a^{n-1}$, to estimate the volume elements and density at the current iteration n:

$$\begin{array}{c}\hfill {\rho }_a^n=\frac{m_a{\sum }_b{X}_b^{n-1}{W}_{\mathit{ab}}(h_a)}{X_a^{n-1}}\mathrm{with}{X}_a^{n-1}={\left(\frac{m_a}{{\rho }_a^{n-1}}\right)}^p.\end{array}$$(14)

The second is to make use of the density calculated in the standard way in the current time-step, ${\rho }_a^0={\sum }_b{m}_b{W}_{ab}(h_a)$, to estimate X_a and then ${\rho }_a^n=m_a/V_a^n$ ($V_a^n$ calculated with Eq. (12)). That is:

$$\begin{array}{c}\hfill X_a={\left(\frac{m_a}{{\rho }_a^0}\right)}^p.\end{array}$$(15)

We refer hereafter to the first method as an explicit procedure driven by X_1, a (Eq. (14)), an to the second as an implicit method, driven by X_2, a (Eq. (15)).

In the simple, static, toy models discussed below, X_1, a convergence is achieved after a few iterations. Once the estimator has converged, there is always an enhancement of the partition of unity, which is almost perfect when the exponent p → 1. Thanks to its simplicity and good results, X_1, a was the estimator chosen in Paper I. It has, however, a couple of drawbacks that are worth noting. Firstly, because of the explicit nature of Eq. (14), a complete Lagrangian consistency is never achieved. Secondly, taking p = 1 makes it too sensitive to particle noise and not recommended. Therefore, p = 0.7 − 0.8 are the recommended values, which slightly degrades the performance of X_1, a and introduces an undesirable free parameter (the particular value of the exponent p). These issues can be overcome with the second procedure to implement X_a which, in the end, was the one chosen in this work. Therefore, X_a consistently refers to X_2, a hereafter.

In the second method, the calculation of the VEs makes use of the value of ${\rho }_a^0$ calculated in the current integration step. As we go on to see in the static toy models (Sect. 4), the explicit option X_1, a leads to a better partition of the unity and interpolations than X_2, a for identical values of the p exponent, but the former is less numerically robust and not fully compatible with the Lagrangian formulation of the SPH equations. On the contrary, the estimator X_2, a allows us to build a Lagrangian-consistent scheme (Appendix A), which incorporates the grad-h terms – provided that the exponent p in Eq. (15) is chosen equal to one (p = 1). This allows us to eliminate a parameter.

In the following sections, we show that reducing the error in the kernel normalization (E₁, hereafter) of particle a, as in:

$$\begin{array}{c}\hfill E_1={\displaystyle [\sum _b{V}_b{W}_{\mathit{ab}}(h_a)-1],}\end{array}$$(16)

usually improves the requirement:

$$\begin{array}{c}\hfill E_2·h_a={|⟨\mathbf{\Delta }\mathit{r}⟩|}_a={\displaystyle \left|\sum _b{V}_b({\mathit{r}}_b-{\mathit{r}}_a)W_{\mathit{ab}}(h_a)\right|\simeq 0,}\end{array}$$(17)

where E₂ is the normalized error module |⟨Δr⟩|_a/h_a of particle a. The error reduction is the consequence of using the estimators X_1, 2 to evaluate the volume elements. Any reduction in both errors (E₁, E₂) will potentially improve the dynamic evolution of the simulated physical system.

4. Estimating the errors E₁ and E₂ for different particle distributions

A good control of errors E₁ and E₂ is of utmost importance to the SPH technique. This is because the quality of both the interpolated function ⟨f(r)⟩,

$$\begin{array}{c}\hfill ⟨f(\mathit{r})⟩\simeq f(\mathit{r}){\displaystyle \sum _b{V}_b{W}_{\mathit{ab}}(h_a)+\mathbf{\nabla }\mathit{f}·\sum _b{V}_b({\mathit{r}}_b-{\mathit{r}}_a)W_{\mathit{ab}}(h_a),}\end{array}$$(18)

and its gradient ⟨∇f(r)⟩ (Eq. 7), are very sensitive to these errors (Rosswog 2015a). Nevertheless, both errors are correlated, as it can be inferred from the following argumentation in one dimension. We first take the kernel normalization condition as a function of the spatial coordinate,

$$\begin{array}{c}\hfill G(x_a)={\displaystyle \sum _b{V}_b{W}_{\mathit{ab}}(h_a).}\end{array}$$(19)

Using the Gaussian kernel $W_{\mathit{ab}}(h_a)=\frac{C}{h_a}exp[-{(\frac{x_b-x_a}{h_a})}^2]$, the standard SPH derivative of G(x) is expressed as:

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}G}{\mathrm{d}x}\right)}_a=-\frac{2}{h_a^2}{\displaystyle \sum _b{V}_b(x_b-x_a)W_{\mathit{ab}}(h_a),}\end{array}$$(20)

thus,

$$\begin{array}{c}\hfill {\displaystyle \sum _b{V}_b(x_a-x_b)W_{\mathit{ab}}=-\frac{h_a^2}{2}{\left(\frac{\mathrm{d}G}{\mathrm{d}x}\right)}_a.}\end{array}$$(21)

We note that the LHS of Eq. (21) is in fact E₂, suggesting that having a good partition of the unity (G ≃ 1, i.e., E₁ → 0) makes dG/dx ≃ 0 and, as a consequence, the E₂ error is suppressed. An independent proof of the link between E₁ and E₂ obtained with an exponential kernel (Fulk & Quinn 1996) is given in the Appendix A of Paper I.

It should be recognized, however, that the proof above is only indicative because, unlike the interpolators widely used in practical calculations, the Gaussian and the exponential kernels are functions without compact support. Moreover, it could happen that even when G(x) is close to one, it shows fluctuations around the particles and its derivative may significantly differ from zero.

In this regard, additional insights about the impact of the VE choice in the errors E₁ and E₂ can be obtained by studying a handful of static particle distributions and using kernels with compact-support. These errors are expected to be large close to discontinuities, which in SPH are usually spread over a few times the smoothing-length distance (h). We chose three representative discontinuities that often appear in practical calculations: a Gaussian (model A), an inverted-Gaussian (model B), and a wall (model C). These are given by the following mathematical expressions:

$$\begin{array}{c}\hfill \rho (x)={\rho }_0+\mathrm{\Delta }\rho e^{-{\left(\frac{x-x_0}{\delta }\right)}^2},(A)\end{array}$$(22)

$$\begin{array}{c}\hfill \rho (x)={\rho }_0-\mathrm{\Delta }\rho e^{-{\left(\frac{x-x_0}{\delta }\right)}^2},(B)\end{array}$$(23)

$$\begin{array}{c}\hfill \rho (x)={\rho }_0+\mathrm{\Delta }\rho \frac{e^{\left(\frac{x-x_0}{\delta }\right)}-e^{-\left(\frac{x-x_0}{\delta }\right)}}{e^{\left(\frac{x-x_0}{\delta }\right)}+e^{-\left(\frac{x-x_0}{\delta }\right)}},(C)\end{array}$$(24)

where the values of the parameters ρ₀, Δρ, and δ are specified in Table 1.

We arranged these profiles into a 1D distribution of N = 100 particles with equal mass and distributed according to the density profile, with reflective boundary conditions. We used an interpolating sinc-kernel with exponent n = 5 (Cabezón et al. 2008), which has a shape similar to the M₆ spline. The different profiles of errors E₁ and E₂ at each point, for the estimators X₁ and X₂, are depicted in Fig. 1. As we can see, X₁ leads to a clear improvement in the kernel normalization condition as p increases. For p ≃ 1 the E₁ error becomes negligible, hence it is not shown. Interestingly, the estimator X₂ with p = 1 (black crosses) leads to similar E₁ and E₂ errors as X₁ with p = 0.8 in all profiles. The normalized error E₂ = |⟨Δx⟩|/h follows a similar trend, suggesting that having a good partition of unity is not only beneficial to approaching an estimation of the density, but also to calculate the gradient of any magnitude of interest with ISPH.

Fig. 1. Results of the SPH evaluation of the test profiles A, B, and C using different VEs. Top row: density profiles of models A (Gaussian), B (Inverted-Gaussian), and C (Wall), of Table 1, calculated with X1 using p = 0.0 (magenta) and p = 0.8 (blue) in Eq. (13). Crosses (×) in black are for X2 (Eq. (15), with p = 1) and points in light blue are the analytic values. Central row: E1 (partition of unity) error for the same three profiles following the same color schemes, with the black lines for X2 with p = 1. Bottom row: same as central row but for error E2 = |⟨Δx⟩|/h.

4.1. Impact of VE in estimating gradients

We can use our simple sharp profiles above to gain insight into the relationship between E₁, E₂ in Eqs. (16), (17), as well as the accuracy of the first derivative. To do so, we assume that the density of the test particle distribution follows profile A, Eq. (22), so that it totally determines the VEs through Eq. (13). Let us also assume that we wish to obtain the SPH derivative of a generic wall-like function f given by profile D in Table 1, Eq. (24)⁴. Such a derivative, $\frac{\mathrm{d}f}{\mathrm{d}x}$, is sensitive to the choice of the estimator X to compute the VEs. We can thus compare the analytic and the numerical value of $\frac{\mathrm{d}f}{\mathrm{d}x}$ and carry out the L₁ analysis of the results with Eq. (26).

Figure 2 depicts the wall-like function f(x) (black solid line) and its analytical gradient of df/dx (light-blue line). Although the value around the maximum of $\frac{\mathrm{d}f}{\mathrm{d}x}$ is similar for p = 0.0 (magenta solid line) and p = 0.9 (blue solid line), the last option fits better the derivative around the coordinate x = 0.6, which is in turn similar to that obtained using the estimator X₂ with p = 1 (black dots).

Fig. 2. Wall-like function f(x) (black line) and its analytical gradient df/dx (light-blue solid line). The solid lines in magenta and blue were obtained with X1 and p = 0.0, 0.9 respectively, while the black points are for X2 with p = 1.

The left panel in Fig. 3 shows the averaged L₁ value for E₁ (Eq. (16)) and E₂ (Eq. (17)):

$$\begin{array}{c}\hfill L_1(E_{1,2})=⟨|E_{1,2}|⟩,\end{array}$$(25)

Fig. 3. L1 calculation of the errors in the numerical experiments shown in Fig. 2. Left: averaged value of L1 for the partition of the unity (solid lines) and ⟨Δx⟩/h (dotted lines). Blue lines are for the estimator X1 and magenta lines for X2. Right: averaged value of L1, obtained with Eq. (26), for the derivative of the wall function calculated with X1 (solid blue line) and X2 (solid magenta line). The black line is the calculation with the non-conservative ISPH, which is independent of the particular value of the exponent p. We note that when p ≃ 1 both schemes (ISPH and ncISPH) converge, albeit faster for estimator X1 than for X2.

calculated in the interval 0.3 ≤ x ≤ 0.7 for the test presented above. The L₁ values for the partition of the unity and ⟨Δx⟩/h decreases as the exponent p increases, as expected. Nevertheless, the quantitative details again depend on the type of estimator, X₁ or X₂. As the figure shows, there is a factor ten reduction of the L₁ errors for X₁ and p ≃ 0.9. The errors become negligible when p ≃ 1. Although the option X₂ shows a lower convergence rate at larger p values, it still provides a significant reduction of L₁, decreasing the errors in almost a factor ten when p = 1.

The panel on the right in Fig. 3 shows the averaged L₁ error of $\frac{\mathrm{d}f}{\mathrm{d}x}$ normalized to $⟨{\left(\frac{\mathrm{d}f}{\mathrm{d}x}\right)}_{\mathrm{analytic}}⟩$,

$$\begin{array}{c}\hfill L_1=\frac{1}{N⟨{\left(\frac{\mathrm{d}f}{\mathrm{d}x}\right)}_{\mathrm{analytic}}⟩}{\displaystyle \sum _{b=1,N}}\left|{({\left(\frac{\mathrm{d}f}{\mathrm{d}x}\right)}_{\mathrm{analytic}}-{\left(\frac{\mathrm{d}f}{\mathrm{d}x}\right)}_{\mathrm{sph}})}_b\right|,\end{array}$$(26)

and how both schemes (ISPH and ncISPH) converge in light of the flavor used for the estimator X (X₁ solid blue line and X₂ solid magenta line). We note that because the ncISPH scheme (Eq. (5)) makes exact linear interpolations, the results (black line) are not sensitive on the adopted value of p. On the other hand, the results of the conservative scheme, ISPH (blue and magenta lines), show a clear dependence on the choice of the volume element, as expected. The profile of L₁(p) decreases linearly, achieving the same accuracy as the ncISPH scheme when p > 0.9 for X₁ or close to it for X₂.

Even though the test cases presented in this section were calculated in 1D with ordered particle settings, the results unambiguously support the idea that: a better partition of unity improves both the estimation of density and the calculation of gradients in ISPH. Hereafter, we focus on the volume elements calculated with the estimator X₂ (Eq. (15)) with p = 1, because, unlike X₁, it fits in with the Lagrangian formulation of the SPH equations perfectly and also allows us to eliminate the parameter, p. Additionally, the grad-h terms can be easily incorporated to the scheme with this choice of VEs. The resulting ISPH equations are described in Sect. 5.1.

4.2. The quest for fully implicit VEs

In principle, the best procedure for calculating the volume elements would be to directly obtain them from the inversion of the kernel normalization matrix. Unlike the indirect methods described in the previous section, a fully implicit implementation of the VEs has the advantage of always leading to the perfect partition of the unity. Previous attempts to build implicit SPH schemes (Knapp 2000; Escartín 2016), have made use of advanced techniques to efficiently invert large sparse matrices, as, for example, the PARDISO library⁵. These types of libraries could also be used to implicitly find the volume elements V_b by solving large linear systems with N × N equations and unknowns:

$$\begin{array}{c}\hfill {\displaystyle \sum _{b=1}^{n_b}{V}_b{W}_{\mathit{ab}}(h_a)=1.}\end{array}$$(27)

Unfortunately, finding the VEs with such a direct approach generally leads to non-physical results. We calculated, in a fully implicit manner, the VEs of the Gaussian function (Eq. (22)) by solving the linear system above and found that the volume elements strongly oscillate around the explicit solution given by Eq. (14) with p = 1 (see Fig. 4). Furthermore, in many points, the volume elements become negative, which is also non-physical. Therefore, we discarded the fully implicit route to finding the VEs in this work.

Fig. 4. Implicit versus explicit estimation of the volume elements, Vb, in the numerical experiments with the Gaussian curve (Eq. (22)).

5. ISPH equations

Here, we summarize the set of ISPH equations used to compute the 3D tests on the following section. We provide two sets of momentum and energy equations, one corresponding to a “standard” of generally used VEs (i.e., X_a = m_a) and another one to our suggested VEs, which are the equations used to perform the comparisons in the tests of Sect. 6. We also provide the details of the procedure used to build these equations in Appendix A.

5.1. Description of the ISPH equations

The density equation is expressed as:

$$\begin{array}{c}\hfill {\rho }_a=\frac{m_a}{V_a},\end{array}$$(28)

with V_a defined by Eq. (12).

The momentum equation is:

$$\begin{array}{c}\hfill {\ddot{x}}_{i,a}=\{\begin{array}{c}\mathrm{For}{X}_{a,b}=m_{a,b}:\hfill \\ -\frac{X_a}{m_a}{\sum }_b{m}_b[\frac{X_b{P}_a}{{\mathrm{\Omega }}_a{k}_a^{2-\sigma }{k}_b^{\sigma }}{\mathcal{A}}_{i,ab}(h_a)+\frac{X_b{P}_b}{{\mathrm{\Omega }}_b{k}_b^{2-\sigma }{k}_a^{\sigma }}{\mathcal{A}}_{i,ab}(h_b)];\hfill \\ \hfill \\ \mathrm{For}{X}_{a,b}=m_{a,b}/{\rho }_{a,b}^0:\hfill \\ -{\sum }_b{m}_b[\frac{X_a^{2-\sigma }{X}_b^{\sigma }{P}_a}{{\mathrm{\Omega }}_a{m}_a^2{k}_a}{\mathcal{A}}_{i,ab}(h_a)+\frac{X_b^{2-\sigma }{X}_a^{\sigma }{P}_b}{{\mathrm{\Omega }}_b{m}_b^2{k}_b}{\mathcal{A}}_{i,ab}(h_b)];\hfill \end{array},\end{array}$$(29)

with ${\rho }_a^0={\sum }_b{m}_b{W}_{ab}(h_a)$, and Ω_a, b is given in Appendix A. We introduced a parameter, 0 ≤ σ ≤ 1, which allows us to choose between the pure Lagrangian scheme developed in the Appendix A (σ = 0) and a progressive deviation of it, which suppresses the tensile instability when σ → 1 (see Sect. 5.2 for more details on σ). Moreover,

$$\begin{array}{c}\hfill {\mathcal{A}}_{i,ab}(h_{a,b})={\displaystyle \sum _{j=1}^d{c}_{ij,a}(h_a)(x_{j,b}-x_{j,a})W_{\mathit{ab}}(h_{a,b}),}\end{array}$$(30)

with c_ij, a as the coefficients of the inverse matrix in the IA given by Eq. (7), and d as the number of dimensions. Any expression of standard SPH can indeed made compatible with the IA by taking the kernel derivative as (Cabezón et al. 2012):

$$\begin{array}{c}\hfill \frac{\partial W_{\mathit{ab}}(h_a)}{\partial x_{i,a}}={\mathcal{A}}_{i,ab}(h_a);i=1,d.\end{array}$$(31)

The energy equation is:

$$\begin{array}{c}\hfill {\dot{u}}_a=\{\begin{array}{c}\mathrm{For}{X}_{a,b}=m_{a,b}:\hfill \\ \frac{X_a{P}_a}{m_a{\mathrm{\Omega }}_a{k}_a^{2-\sigma }}{\sum }_b{\sum }_{i=1}^d\frac{X_b}{k_b^{\sigma }}\left[(v_{i,a}-v_{i,b}){\mathcal{A}}_{i,ab}(h_a)\right];\hfill \\ \hfill \\ \mathrm{For}{X}_{a,b}=m_{a,b}/{\rho }_{a,b}^0:\hfill \\ \frac{X_a^{2-\sigma }{P}_a}{m_a^2{\mathrm{\Omega }}_a{k}_a}{\sum }_b{\sum }_{i=1}^d{m}_b{X}_b^{\sigma }\left[(v_{i,a}-v_{i,b}){\mathcal{A}}_{i,ab}(h_a)\right];\hfill \end{array}.\end{array}$$(32)

The necessary terms of artificial viscosity (AV) are added to the right of the equations above as in Paper I, but explicitly including a quadratic term in the signal velocity:

$$\begin{array}{c}\hfill v_{\mathit{ab}}^{\mathrm{sig}}={\overline{\alpha }}_{\mathit{ab}}{\overline{c}}_{s,ab}-\beta w_{\mathit{ab}},\end{array}$$(33)

where ${\mathit{w}}_{\mathit{ab}}={\mathit{v}}_{\mathit{ab}}·{\widehat{\mathit{r}}}_{\mathit{ab}}$ (Monaghan 1997; Price et al. 2018). The parameter β is kept constant with a default value β = 2. The AV coefficient α is controlled with the switches scheme described in Read et al. (2010) so that α ≃ 1 in strong shocks but it decays to a minimum value of α ≃ 0.05 away from them. Finally, ${\overline{\alpha }}_{\mathit{ab}}$ and ${\overline{c}}_{s,ab}$ are the average α and speed of sound between neighboring particles. According to Monaghan (1997), the AV contribution to the energy equation should include a heat-conduction term which smooths the pressure in wall-shock conditions. The precise form of such heat-conduction term is:

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}{u}_a}{\mathrm{d}t}\right)}_{\mathrm{cond}}^{\mathit{AV}}={\displaystyle \sum _b\sum _{i=1}^d{m}_b{\alpha }_u\frac{v_{ab,\mathrm{cond}}^{\mathrm{sig}}(u_a-u_b)}{{\overline{\rho }}_{\mathit{ab}}}\frac{r_{i,ab}}{|r_{\mathit{ab}}|}{\overline{\mathcal{A}}}_{i,ab}(h_a,h_b),}\end{array}$$(34)

with $\overline{\mathcal{A}}(h_a,h_b)=0.5(\mathcal{A}(h_a)+\mathcal{A}(h_b))$, which in the SPHYNX code is implemented as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{\mathrm{d}{u}_a}{\mathrm{d}t}\right)}_{\mathrm{cond}}^{\mathit{AV}}=& {\displaystyle \sum _b\sum _{i=1}^d\frac{1}{2}{\alpha }_u{v}_{ab,\mathrm{cond}}^{\mathrm{sig}}(u_a-u_b)}\hfill \\ \hfill & \{V_a\frac{m_b}{m_a}\frac{r_{i,ab}}{|r_{\mathit{ab}}|}{\mathcal{A}}_{i,ab}(h_a)+V_b\frac{r_{i,ab}}{|r_{\mathit{ab}}|}{\mathcal{A}}_{i,ab}(h_b)\}\hfill \end{array}.\end{array}$$(35)

The signal velocity ${\mathit{v}}_{\mathrm{cond}}^{\mathrm{sig}}$ used in our tests is (Price 2008):

$$\begin{array}{c}\hfill v_{ab,\mathrm{cond}}^{\mathrm{sig}}=\sqrt{\frac{|P_a-P_b|}{{\overline{\rho }}_{\mathit{ab}}}}.\end{array}$$(36)

The results of the tests below suggest that adding a small amount of conductive heat by Eq. (35) is beneficial because it contributes to reduce the tensile instability and to smooth the numerical noise. Nonetheless, the value of the constant α_u should not be high, otherwise the density-peak in strong shocks may be under estimated (see Sect. 6.4). We chose α_u = 0.1, which is low and in agreement with the choice by other authors, such as Tricco (2019). The equation of state (EOS) is that of an ideal gas, with γ = 5/3.

5.2. Obtaining the value of σ

When using the SPH equations deduced from the Euler–Lagrange formulation a drawback appears wherever ∇ρ becomes large as in the case of contact discontinuities, for instance. In such cases, the incorrect estimation of gradients may lead to numerical artifacts, with the tensile instability standing as one of the most harmful and most common. Several solutions have been postulated to cope with this problem, all of them requiring some departure from the exact Lagrangian formulation. For example, Ritchie & Thomas (2001) proposed to estimate the density by averaging over the internal energies so that the ensuing density field is smooth. A similar approach was described in Saitoh & Makino (2013), who suggested that the volume elements be redefined so that they depend on the pressure rather than density. Another solution was proposed by Read et al. (2010), where a typical element within the summations of the momentum and energy equations is changed from the standard $[P_a/{\rho }_a^2]$ to [P_a/(ρ_aρ_b)]. While it is not totally Lagrangian compatible, such a simple change is mathematically consistent with the standard derivation of the SPH equations (Monaghan 1992; Price 2004)⁶ and can totally suppress the tensile instability. All of these SPH variants have a feature in common, namely, that the main magnitudes in the movement equation are somehow “crossed” for particles a and b (f.e. u_b/ρ_a in Ritchie & Thomas 2001 or P_a/(ρ_aρ_b) in Read et al. 2010). In Cabezón et al. (2017), the choice of X_1, a (Eq. (14)), combined to the SPH equations defined in Rosswog (2015b), naturally led to an almost “crossed” scheme⁷. In spite of not being totally Lagrangian-compatible and with the grad-h terms only approximated, this code behaved well and was able to pass a suite of verifying tests.

Here, we propose a similar procedure as in Read et al. (2010), but allowing for a self-adaptive crossing of indexes in the original Lagrangian equations derived in the Appendix A. The resulting expressions, Eqs. (29) and (32), are steered by a parameter σ so that 0 ≤ σ ≤ 1. The value σ = 0 leads to the fully Lagrangian SPH equations whereas σ = 1 gives a fully crossed expression. We make use of a ramp function (R) to automatically decide the instantaneous value of σ as a function of the density contrast between any pair of particles.

$$\begin{array}{c}\hfill {\sigma }_{\mathit{ab}}=\{\begin{array}{cc}0;\hfill & A{t}_{\mathit{ab}}\le A{t}_{min}\hfill \\ R(A{t}_{\mathit{ab}}-A{t}_{min});\hfill & A{t}_{min}\le A{t}_{\mathit{ab}}\le A{t}_{max}\hfill \\ 1;\hfill & A{t}_{\mathit{ab}}\ge A{t}_{max}\hfill \end{array},\end{array}$$(37)

with $R=\frac{1}{A{t}_{max}-A{t}_{min}}$ and $A{t}_{\mathit{ab}}=\left|\frac{{\rho }_a-{\rho }_b}{{\rho }_a+{\rho }_b}\right|$ as the Atwood number. Equation (37) leads to σ ≃ 1 only in those regions hosting large density gradients, while the fully Lagrangian scheme is preserved wherever σ = 0, which is taken by the majority of the particles of the system. We empirically found that At_min = 0.1, At_max = 0.2 produce satisfactory results in all the numerical experiments described in this work.

6. Tests

In this section, we analyze the results of four well-known tests that require the solution of the full system of Euler equations in three dimensions. The first test deals with the hydrostatic equilibrium of a two-phase fluid. It aims to analyze the abilities of our Lagrangian scheme to handle sharp density discontinuities. Then we applied this scheme to the study of the growth of the Kelvin–Helmholtz (KH) instability, to the interaction of a supersonic wind with an isothermal, high-density spherical cloud (the wind-cloud test), as well as to simulate the evolution of a point-like explosion (the Sedov test). The simulations were carried out using the ISPH code SPHYNX and the outcomes have been compared to well-known solutions. Tests were run with varying parameters to disentangle the effects of IA, VE, heat transport in the AV, and σ. Special emphasis is placed on a comparison among models with different choices of σ and volume elements, namely, the standard choice⁸ of X_a = m_a and the improved VEs with $X_a=m_a/{\rho }_a^0$.

We find that the conservation properties are excellent and that linear momentum is preserved to almost machine precision in all tests. The angular momentum, with respect to the center of mass, deviates less than 10⁻⁶ from the initial value. The relative deviation, |ΔE|/E₀ in the total energy was less than 10⁻⁵, except in the isobaric two-fluid test (Sect. 6.1), where it was ≃10⁻⁴ at late times.

6.1. Isobaric two-fluid test

The simulation of the hydrostatic evolution of a two-phase system with very different densities and internal energies is far from being trivial with SPH codes. We considered a system of two fluids separated by a contact discontinuity but in pressure equilibrium:

$$\begin{array}{c}\hfill \rho =\{\begin{array}{cc}4\hfill & 0.25\le x,y,z\le 0.75,\hfill \\ 1\hfill & \mathrm{otherwise}.\hfill \end{array}\end{array}$$(38)

The system is isobaric with P = 2.5 and we use N = 110³ equal-mass particles spread in two nested body-centred cubic (bcc) lattices. The EOS is P = (γ − 1)ρu with γ = 5/3, the number of neighbors is set to n_b = 100, and the sinc (n = 5) kernel (Cabezón et al. 2008) is used to interpolate. If the density around the contact discontinuity is not adequately smoothed, the two phase system evolves non-physically when the full Lagrangian scheme (Eqs. (29) and (32), with σ = 0) is applied. The reason is that the error in ∇ρ becomes too large at the contact discontinuity, inducing the tensile instability. This is just the kind of simulation where the use of the magnitude σ (see Sect. 5.2), defined by Eq. (37), becomes especially helpful.

The results of applying Eqs. (29) and (32) to the hydrostatic square test are summarized in Table 2 and Fig. 5. Only models with $X_a=m_a/{\rho }_a^0$ are shown in Fig. 5, and calculations with X_a = m_a give similar results, albeit with higher L1 values. The density color map slices depicted in Fig. 5 show that the behavior of the square is primarily controlled by the value of σ. In the case of the full Lagrangian formulation (models H₅, H₆), calculated with σ = 0, the system completely looses its shape in half of the sound-crossing time, τ_sc ≃ 0.9 s. Nevertheless, increasing the conductivity in the AV (model H₆) slightly delays the deformation; although, in the end, it is not able to prevent it. Model H₇, calculated with σ = 1, leads to the best results. The cube maintains its profile at t/τ_sc ≃ 1.7 (last snapshot in third row in Fig. 5) and beyond. Such behavior is in fair agreement with the results of other density-based schemes (Read et al. 2010; Wadsley et al. 2017) calculated in 2D and 3D respectively. It also matches the 2D calculation in Paper I (first row in Fig. 4 of that work, with 0 ≤ t/τ_sc ≤ 2) calculated with the explicit X_1, a given by Eq. (14) with p = 0.7 (see footnote 7). It is worth noting that model H₈, calculated with σ[0 − 1], preserves the shape of the square until t/τ_sc ≃ 1 and still does at the last snapshot in Fig. 5 at t/τ_sc = 1.7, although not as well as in model H₇. As shown in Fig. 6 (upper panel), the algorithm (Eq. (37)) that self-adapts σ as a function of the local density contrast works splendidly, increasing the value of σ only where it is needed.

Fig. 5. Slices around Z = 0 showing the density colormap of models H5, H6, H7, H8 (Very similar results for H1–H3 models, calculated with standard VEs can be found in plots given in previous works (Cabezón et al. 2017, their Fig. 4).) in Table 2 (rows) at different times (columns), with the sound-crossing time being τsc = 0.9.

Fig. 6. Isobaric two-fluid numerical experiment. Upper panel: slice around Z = 0 showing the averaged σ-parameter (Eq. (37)), in model H8 at t = 0.55. Appreciable values of σ are only attained at the fluid inter-phase. Lower panel: time evolution of averaged L1(E1) (solid lines), L1(E2) (dashed lines) for models H3, H7, and H8 in Table 2.

To investigate the dependence of the errors E_1, 2 (Eqs. (16) and (17)) with respect to the VE estimator and σ we calculated the L₁ values (Eq. (25)) in a rectangular shell with width 0.3 around the contact discontinuity.

The column L₁(d) of Table 2 is defined by the average error of the absolute displacement of the SPH particles located in the shell, with respect to their initial positions:

$$\begin{array}{c}\hfill L_1^d(t)=\frac{1}{N}{\displaystyle \sum _{b=1}^N\sqrt{{(x_b(t)-x_{b,0})}^2+{(y_b(t)-y_{b,0})}^2+(z_b(t)-z_{b,0}^2)},}\end{array}$$(39)

where x, y, and z are the Cartesian coordinates of the particles while x₀, y₀, and z₀ coordinates represent their initial position. Sub-index b runs from 1 up to the N particles contained in the rectangular shell.

The full Lagrangian models H₁, H₂, H₅, and H₆, calculated with σ = 0, give the worst results, with L₁(d)≥3 × 10⁻² at t = 1.5. The minimum value of L₁(d) corresponds to model H₇ (L₁(d) = 1.5 × 10⁻²), calculated with σ = 1 and improved VEs. However, model H₈ calculated with the self-adaptive σ[0 − 1] and improved VEs, also displays a similarly good behavior, while keeping full Lagrangian consistency (σ ≃ 0) in a large fraction of the domain (see Fig. 6, upper panel).

The lower panel in Fig. 6 depicts the evolution of the averaged estimator L₁ of the partition of the unity (solid lines) and ⟨Δr⟩/h (dashed lines) for models H₃, H₇, and H₈ in Table 2. Comparing H₃ to H₇ highlights the relevance of the VEs, while comparing H₇ to H₈ shows the impact of a variable σ. Is it clear that the partition of the unity is enhanced with the choice X = m/ρ⁰, nevertheless the behavior of ⟨Δr⟩/h is less defined. The error L₁(E₂) remains lower for models H₇ and H₈ than model H₃ until t ≃ 0.2. Afterwards the error grows to eventually reach a similar value to that of model H₃, calculated with X_a = m_a. Thus, we expect a degradation in the behavior ⟨Δr⟩/h with respect the results shown in Figs. 1 and 3 whenever the disorder of the particles is high. In summary, the dominant parameter for the dynamic evolution of the system is σ, which basically determines the stability of the square. Using X = m/ρ⁰ improves the partition of unity at any time, but the magnitude ⟨Δr⟩/h does not show a definite enhancement in this test.

6.2. Kelvin–Helmholtz instability

The correct simulation of the evolution of the contact layer between fluids with different densities is of capital importance for the adequate growth of the Kelvin–Helmholtz (KH) instability. García-Senz et al. (2012) proved that the use of ISPH improves the evaluation of gradients overall (particularly in the contact layer). Provided the density contrast is not very high, it prevents the appearance of tensile instabilities that otherwise suppress the growth of the KH instability. Because SPHYNX uses ISPH by default, all KH simulations presented here show no signs of tensile instability and have growth rates close to the reference calculation by McNally et al. (2012) with the code PENCIL. Nevertheless, using a volume element that is better at fulfilling conditions (16) and (17) should additionally improve the accuracy of ISPH and, as a consequence, obtain a better KH growth rate.

We ran this test in a thin three-dimensional layer of size [1 × 1 × 0.0625] with 4.2 × 10⁶ equal-mass particles, distributed in a ratio 2:1 between the high and low density regions. For the initial setting, we have three stratified layers, being the central layer the high-density one. Each region was generated from a random particle distribution relaxed to a glass-like configuration. The equation of state, initial velocities, and initial pressure are the same as those in Sect. 5.4.1 in Paper I. Times are normalized to the characteristic Kelvin–Helmholtz growth time t_KH, as defined in Agertz et al. (2007), which results in t_KH = 1.06 for our models. As discussed at the end of Sect. 6.4, we adopted a low value, namely, α_u = 0.1, of the AV conductive parameter.

We used two different VEs (the standard X_a = m_a and the enhanced version proposed here $X_a=m_a/{\rho }_a^0$) and, for each case, we performed two simulations: one fixing σ = 0, hence ensuring fully Lagrangian compatibility, and another one allowing σ to vary according to Eq. (37). All simulations are summarized in Table 3.

In Fig. 7, we represent the particle distributions for each simulated model (each snapshot corresponds to one model, from KH₁ to KH₄) at t = 2. The color represents density. As it can be seen, there are very few differences among all the simulated models. In all cases, the KH billows are able to develop, with small differences among simulations, mostly at the extremes of the billows.

Fig. 7. Particles distribution in a thin slice around Z = 0 for models KH1 to KH4 at t = 2 (tKH = 1.06). The density [1:2] is color-coded.

We can also explore the growth of the amplitude mode for the v_y field and compare it with a reference evolution taken from the PENCIL code (McNally et al. 2012). Figure 8 (left) shows this evolution and it is clear that there are very little differences among both VEs implementations, in agreement with the results in Paper I. These results are in good match with those by Frontiere et al. (2017) and Rosswog (2020a) for similar resolution and initial conditions, provided that the Cullen & Dehnen (2010) AV trigger is chosen to steer dissipation in the last case.

Fig. 8. Amplitude growth of the vy field in the KH instability test for all calculated models. Solid lines correspond to αmin = 0.05, while dashed lines are from the models with increased AV (αmin = 0.5). Those with a density jump by a factor of 2 are compared with the reference PENCIL simulation in the left panel. Right panel: two simulations with a density jump by a factor of 8 for both choices of the VEs.

Figure 9 presents the time evolution of the averaged L₁ value for errors E₁ and E₂. To average the L₁ values we restricted to those particles that had σ ≥ 0.025. This condition returns all those particles that are found in the inter-phase (see Fig. 10), which is the region the error comes from and where the accurate evaluation of the VEs is more critical. It is clear that whenever the improved VEs are used (dashed lines in Fig. 9), these errors decrease, independently of the scenario being simulated. The last two columns of Table 3 show the numerical value of L₁ for all simulated tests at t = 2, where we can see that the choice X = m/ρ⁰ reduces both, the L₁(E₁) and the L₁(E₂) errors.

Fig. 9. L1 evolution of errors E1 and E2 for the KH models. Solid lines correspond to X = m, while dashed lines are those of X = m/ρ0.

Fig. 10. Color map of σ (Eq. (37)), in a thin slice around Z = 0 for model KH4 at t = 2.

From these results, we can conclude that despite decreasing the errors E₁ and E₂, the VEs choice has little influence in the growth of the KH instability when the density contrast is moderate. The most important element here is having the IA implementation to calculate gradients, as shown in Paper I. Additional improvement in the conditions (16) and (17) via using improved VEs is subdominant. The dominant source of error at this level is the artificial viscosity. To test this we performed two additional simulations, KH₅ and KH₆, with an artificially increased α_min in our AV switches. This α_min controls the amount of dissipation when the AV switch is off. Therefore, is the minimum dissipation present in the whole system. Our standard value is α_min = 0.05, while for KH₅ and KH₆ we used α_min = 0.5. To have a reference, traditional AV without switches employ a global value α_min ≃ 1 − 1.5. The amplitude growth of the v_y field is also shown in the left panel of Fig. 8 (dashed lines). It is clear that the AV has a major effect in the development of the instability. Still, the VEs have a negligible impact when the IA formalism is used.

Our final test was to explore if the VEs can exert an influence when the density contrast is higher. To do that we decreased the density of the low-density regions by a factor of 4, simply by generating a new relaxed model with four times fewer particles, along with repeated simulations KH₂ and KH₄, now with this new density jump by a factor of 8 (simulations KH₇ and KH₈). We present the evolution of the amplitude growth of the v_y field in the right panel of Fig. 8. In this case, the growth is more irregular and there is a clear difference between the VEs in the linear regime. The improved VEs are able to growth faster than the standard VEs in the initial stages of the development of the KH instability, despite converging later, in the non-linear phase. The IA formalism still does a good job evaluating gradients even with a density jump by a factor 8. Nevertheless, using improved VEs at this density contrast shows to have a noticeable effect.

In summary, if there is a mild density contrast, the IA implementation is good enough to make the choice of VEs sub-dominant. The most relevant parameter in these situations is to keep the dissipation at its lowest possible value, so that random noise is still dissipated, but it doesn’t affect the shear between the different fluid layers. In this respect, an improved handling of the AV would be welcomed, such as that presented in Rosswog (2020a) based on the instantaneous numerical entropy generation rate. If there is a bigger density contrast, the IA formalism can be further improved by using a VE that better fulfills conditions (16) and (17) along with σ[0 − 1].

6.3. Wind-cloud collision

The wind-cloud collision scenario, also called “blob” test (Agertz et al. 2007) is a challenging test for SPH, involving several pieces of physics such as strong shocks and mixing due to the KH instability, amidst a two-phase medium with a large density contrast. In this test, a spherical cloud of cold gas, initially at rest, is swept by a low-density stream of gas (the wind) moving supersonically. As a consequence, the cloud deforms and, after a while, it is fragmented and mixed with the background owing to the combined effect of ablation and hydrodynamic instabilities. This scenario has been extensively studied during the last years (Agertz et al. 2007; Read et al. 2010; Saitoh & Makino 2013; Hopkins 2013; Frontiere et al. 2017), where the tensile instability appearing at the wind-cloud contact was identified as the main difficulty to overcome.

We sought to check the ability of our numerical scheme to handle this test. The initial conditions were those described in Paper I, except that the calculation is now carried out in 3D with approximately N_w = 11.23 × 10⁶ and N_c = 1.23 × 10⁵ particles for the wind and the cloud, respectively. The initial particle distribution for both was that of a stable glass-like configuration. The box has a size {0.25 × 0.25 × 1} with periodic boundary conditions. The high-density cloud is initially located at (0.125, 0.125, 0.125) with a radius R = 1/40 and a density ρ_c = 10, ten times bigger than the surrounding wind. The internal energy of the wind and the cloud are u_w = 3/2 and u_c = 3/20, respectively, so that both regions are in pressure equilibrium. The wind has an initial velocity (2.7, 0, 0). The characteristic KH growth time used in this test to normalize the time is τ_KH = 0.0937.

In Table 4 we present the L₁ values of errors E₁ and E₂ for the simulated models. Both errors are clearly lower when the improved VEs are used, independently of the value of σ. This trend is constant during the whole evolution of the simulation as it can be seen in Fig. 11. In particular, it is clear that the models with σ = 0 lead to the largest L₁ errors in the long run, proving that the pure Lagrangian scheme is the most sensitive to tensile instability with such large density gradients. As in the case of the KH test, we make use of a variable σ to track the L₁ error generated by the particles in the inter-phase. In Fig. 12 we show a series of snapshots where the color represents the value of σ and proves the capabilities of our algorithm to track density jumps.

Fig. 11. L1 values for E1 (solid lines) and E2 (dashed lines) in the wind-cloud collision tests. Using the improved VEs reduces both L1 values for the whole simulation.

Fig. 12. Particle distribution of the wind-cloud test in a thin slice around Y = 0 for model W6 at t/tKH = 0, 1, 2, and 3. The parameter σ is color-coded.

The density color map at t/t_KH = 1.5 is shown in Fig. 13 for both choices of VEs. We see here that, independently of the VE choice, mixing is taking place, which will eventually destroy the cloud – albeit more efficiently if σ ≠ 0. Figure 14 shows the maximum density attained in the cloud during the first stages of the interaction. Because of the strong impact of the wind, a shock-wave is born, which moves through the dense cloud and compresses its matter. The density jumps from the original value of ρ₀ = 10 and tends toward the limiting value ρ = 4ρ₀, which is characteristic of a strong shock in a gas with γ = 5/3. All models show a similar behavior but the profile from models calculated with X = m/ρ⁰, give in general larger values than models calculated with X = m.

Fig. 13. Particle distribution of the wind-cloud test in a thin slice around Y = 0 for models W3 and W6 at tKH = 1.5. Density is color-coded.

Fig. 14. Evolution of the maximum density in the cloud at t/tKH ≤ 1.6 for all models in Table 4.

In order to characterize the destruction of the high-density cloud, we tracked the percent of surviving cloud mass (criteria: ρ ≥ 0.64ρ_c and u ≤ 0.9u_w) in function of t/t_KH in Fig. 15. From this result, we can see that the choice of the VEs is less relevant than using a crossed formulation of the equations (i.e., variable σ or σ = 1). Indeed, the cloud is mostly destroyed in all simulations, and this is mainly due to the integral formulation of the SPH equations combined with the addition of the heat transfer term in the AV equation. Nevertheless, the surviving fraction of the cloud decreases faster and to lower values – provided we do not use σ = 0, independently of the VE choice. The slight delay in the evolution of the surviving fraction with X = m/ρ⁰ simply reflects the fact that the ensuing VEs better track the density peak at the forward shock moving through the cloud. Such a density enhancement results in the slight delay in crossing the above threshold condition ρ ≥ 0.64ρ_c. The disintegration rate of the cloud depicted in Fig. 15 when σ ≠ 0 is in close agreement with that obtained with the CRKSPH scheme by Frontiere et al. (2017).

Fig. 15. Percent of surviving cloud in function of t/tKH for the wind-cloud collision tests.

In summary, when large density contrasts are present, better accuracy is achieved by the proposed generalized volume elements, while the suppression of the tensile instability is handled via a surgical departure of the Lagrangian scheme, using a variable σ.

6.4. Sedov explosion

We simulated this test in a three-dimensional square box of side L = 1 and total number of particles N = 100³ arranged in a glass-like configuration. The mass of all the particles was the same, resulting in an homogeneous initial density profile with ρ = 1 ± 0.005. The explosion was initiated at the center of the box by depositing an amount of energy ΔU = 1. That energy was spread following a Gaussian profile with characteristic width δ = 0.1.

To investigate the dependence of the errors E_1, 2 (Eqs. (16) and (17)) with respect the VEs estimator choice and σ, we calculated the L₁ errors in the shocked region⁹ with Eq. (25).

The outcome of the simulations is summarized in Table 5 and Figs. 16–18. The best results were obtained in those models which incorporate the improved VEs. Unlike in the hydrostatic square test, the particular value adopted for the magnitude σ plays a secondary role in the Sedov test. Figure 16 shows the profiles of density at t = 0.09 for models S₄ and S₈. In both cases, the self-adaptive algorithm to estimate σ was active. The best match with the analytic profile (black line) was obtained with $X_a=m_a/{\rho }_a^0$, not only in the peak values, which were better reproduced, but also in the width of the shock front. The relative errors of the maximum density with respect to the analytical value at t = 0.09 are ≃16% and ≃7% for the standard (S₄) and improved VEs (S₈), respectively. Figure 17 shows the distribution of the averaged σ in a thin slice along the system for model S₈. As it can be seen, the value of σ self-adapts so that it approaches to 1 across the shock front, goes down to ≃0.5 in the post-shock region and vanishes in the central region.

Fig. 16. Sedov explosion for models S4 and S8 in Table 5. Top row: density profiles at t = 0.09 calculated with Xa = ma (model S4, standard VEs) and $X_a=m_a/{\rho }_a^0$ (model S8, improved VEs). The black line is the analytical solution. Note that the spherical symmetry is well preserved in both cases and that the density peak is very well reproduced with the improved VEs. Bottom row: same but with the partition of unit. All particles are represented in the plots.

Fig. 17. Color map slice around Z = 0 showing the value of the averaged σ-parameter defined in Eq. (37), for model S8 in Table 5 at t = 0.09.

Fig. 18. Sedov explosion. Evolution of the averaged L1 error, (Eq. (25)), in the shocked region corresponding to models S4,  (Xa = ma), and S8,  ($X_a=m_a/{\rho }_a^0$), in Table 5. We show both the partition of unity and the normalized ⟨Δr⟩/h condition, for the estimator choices: Xa = ma (purple and orange lines) and $X_a=m_a/{\rho }_a^0$ (blue and magenta lines). Also shown is the evolution of the maximum density for the two volume estimators (dashed lines, black is for S4 and light-blue for S8) during the explosion.

The profile of the normalization function ∑_bV_bW_ab(h_a) at t = 0.09 is also shown in Fig. 16, bottom row. As it can be seen, the choice $X_a=m_a/{\rho }_a^0$ substantially improves the partition of the unity with respect the standard choice X_a = m_a. Such assessment is quantitatively confirmed by the temporal evolution of the error estimators L₁(E₁, E₂) shown in Fig. 18. The enhancement is especially good in the case of the partition of unity, although less pronounced in the ⟨Δr⟩/h condition. Therefore, we conclude that combining ISPH with the VEs obtained with $X_a=m_a/{\rho }_a^0$ reduces both errors and improves the simulations.

The impact of increasing the conductive term in the AV was also analyzed. Raising the parameter α_u from our default choice, α_u = 0.1 to α_u = 0.5, substantially reduces the density peak from ρ_max(t = 0.09) = 3.34 in model S₁ to ρ_max(t = 0.09) = 3.16 in model S₂. In light of these results, it is advisable to choose a not too high value for α_u.

6.5. Choice of σ and entropy evolution in a shock

In light of the results of the tests presented here, the question of why we would need a variable σ at all might arise. Obviously, the case σ = 0 would be the most desirable one, as it is fully Lagrangian-compatible. Nevertheless, it is clear that the crossed σ = 1 scheme leads to a more efficient removal of the tensile instability, but the price to pay is the lack of a complete compatibility between the Euler–Lagrange formulation and the standard density equation (Eq. (28)). Such a weakness with regard to the crossed scheme is highlighted in the following test, which considers the behavior of entropy in weak and strong shocks. We want to prove that the entropy does not evolve as well when σ ≠ 0. We additionally show that the new volume elements have a positive impact on the evolution of entropy.

Our numerical experiment is the same as in the Sedov test above but with an additional case where the initial energy of the explosion has been reduced in a factor of ten to simulate the case of a weak explosion. During the evolution of the explosion, the contribution to the entropy of the gas is separated in two components: that arising from pressure forces Δs_P and that from the artificial viscosity Δs_AV. Ideally, Δs_P ≃ 0 because the dissipation is basically driven by the AV.

The second law of thermodynamics states:

$$\begin{array}{c}\hfill T\mathrm{d}s=\mathrm{d}u-\frac{P}{{\rho }^2}\mathrm{d}\rho ,\end{array}$$(40)

where T, s, u, P are temperature, specific values of entropy and internal energy, and pressure. For an ideal gas with specific heat at constant volume, c_v, we have:

$$\begin{array}{c}\hfill \frac{\mathrm{d}s}{c_{\mathrm{v}}}=\frac{\mathrm{d}u}{u}-(\gamma -1)\frac{\mathrm{d}\rho }{\rho },\end{array}$$(41)

with u = c_vT. Thus, integrating on both sides:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }s}{c_{\mathrm{v}}}=ln\left(\frac{u}{u_0}\right)-(\gamma -1)ln\left(\frac{\rho }{{\rho }_0}\right),\end{array}$$(42)

where ρ₀, u₀ are the initial values of density and internal energy. To obtain the evolution of entropy of the ideal gas, without the AV contribution, we first calculate the increment of entropy (Δs)_AV/c_v due to the AV. Subtracting (Δs)_AV/c_v from (42) gives the ideal gas contribution, which should remain negligible during an adiabatic evolution.

The AV contribution to the entropy in a process involving an energy/heat variation dQ is:

$$\begin{array}{c}\hfill \frac{\mathrm{d}Q}{T}={(\mathrm{d}s)}_{\mathrm{AV}}\to \frac{c_{\mathrm{v}}\mathrm{d}Q}{u}={(\mathrm{d}s)}_{\mathrm{AV}}\to \frac{\mathrm{d}{s}_{\mathrm{Av}}}{c_{\mathrm{v}}}=\frac{\mathrm{d}Q}{u},\end{array}$$(43)

which, integrating on both sides, leads to:

$$\begin{array}{c}\hfill \frac{{(\mathrm{\Delta }s)}_{\mathrm{AV}}}{c_{\mathrm{v}}}={\displaystyle \int }\frac{\mathrm{d}Q}{u}={\displaystyle {\int }_0^t}\frac{1}{u}{\left(\frac{\mathrm{d}u}{\mathrm{d}t}\right)}_{\mathrm{AV}}\mathrm{d}t,\end{array}$$(44)

where ${\left(\frac{\mathrm{d}u}{\mathrm{d}t}\right)}_{\mathrm{AV}}$ is directly obtained from the SPH code. The variation of the entropy of the gas, due to pressure forces (Δs)_P/c_v, and excluding the AV contribution, is determined subtracting Eq. (44) from Eq. (42):

$$\begin{array}{c}\hfill \frac{{(\mathrm{\Delta }s)}_{\mathrm{P}}}{c_{\mathrm{v}}}=\frac{\mathrm{\Delta }s}{c_{\mathrm{v}}}-\frac{{(\mathrm{\Delta }s)}_{\mathrm{AV}}}{c_{\mathrm{v}}}.\end{array}$$(45)

In absence of heat transport and providing that there are neither entropy sources (e.g., nuclear reactions) nor entropy sinks, (Δs)_P/c_v ≃ 0.

In Fig. 19, we show the radial profile of |Δs_P|/c_v for the different calculated models. First of all, we note that there is a spurious variation of entropy in all cases. Nevertheless, the combination $(X_a=m_a/{\rho }_a^0;\sigma =0)$ shown in the fourth panel leads to the best behavior with regard to the entropy in the shocked region (magenta lines). Not surprisingly, the worst cases are for σ = 1 (blue lines). These results clearly prove that Lagrangian compatibility is of utmost relevance for keeping the entropy variations at low levels in adiabatic flows. Furthermore, upgrading the partition of unity with a Lagrangian-compatible estimator, such as $X_a=m_a/{\rho }_a^0$, substantially improves the results. The mixed option with σ = [0 − 1] (black lines) is, for most of the cases, closer to σ = 0 or in between both curves, as expected. It therefore represents a balanced solution between Lagrangian-compatibility and the suppression of the tensile instability.

Fig. 19. Radial profile of the entropy variation by pressure forces, |ΔsP|/cv, for weak (upper row, at t = 0.18) and strong (lower row, at t = 0.08) point-like explosions and different σ choices. The first column is for estimator Xa = ma and the second for $X_a=m_a/{\rho }_a^0$.

The profiles of the total entropy variation is shown in Fig. 20. The strong shock dissipates much more entropy than the weak. Thus, the spurious generation of entropy in weak shocks, as shown in the upper panels of Fig. 19, is comparatively more acute.

Fig. 20. Radial distribution function of total entropy, |Δstot|/cv, released in weak or strong shocks with σ = 0 and the two VEs considered in this work. The entropy in the weak shock cases has been multiplied by 10.

The results of this test support the idea that working with an adaptive σ is a useful option that allows us to suppress the tensile instability across shocks and contact discontinuities, while reducing the unwanted growth of entropy in the post-shock regions.

7. Conclusion

The main goal of this work is to devise and check an SPH scheme that fulfils as many basic physical and numerical constraints as possible, with regard to the demands of modern hydrodynamic codes. The former include conservation of mass, linear and angular momentum, energy, and entropy. The latter requires a good depiction of gradients and shocks besides a correct numerical representation of basic kernel properties, as for example, the kernel normalization. The proposed scheme strongly relies in the ISPH method plus enhanced volume elements described in Cabezón et al. (2017), but with an improved Lagrangian compatibility and a better treatment of sharp density gradients.

To this end, here we introduce an easy scheme to improve the partition of the unity in SPH, with a special emphasis on the connections between the partition of unity and the accuracy in estimating the gradients. By combining analytical reasoning with simple 1D toy models and full 3D simulations, we have shown that improving the constraint ∑_bV_bW_ab(h_a) = 1 automatically leads to an enhancement of the condition ⟨Δr⟩=∑_bV_b(r_b − r_a)W_ab(h_a) = 0. When gradients are approached with an integral expression (the ISPH scheme via Eq. (7)), the fulfillment of the ⟨Δr⟩ constraint is shown to be a sufficient condition to perform exact linear interpolations.

One of the novelties of this work, is the improvement of the partition of the unity, without leaving the Lagrangian formulation of ISPH, so that the enhancement is totally compatible with the inclusion of grad-h terms in the momentum and energy equations. To do that the volume elements, V_a, have been re-defined so that V_a = X_a/k_a with $X_a=m_a/{\rho }_a^0$; $k_a=\sum _b{X}_b{W}_{ab}(h_a)$ and ${\rho }_a^0=\sum _b{m}_b{W}_{ab}(h_a)$. That particular choice of V_a results in the new set of ISPH equations described in Appendix A. These ISPH equations are not only Lagrangian-compatible, displaying perfect conservation properties, but they manifestly improve the partition of the unity which translates into a better estimation of gradients. The computational cost of including these improved VEs, mainly due to the implicit search of the coupled h_a and ${\rho }_a^0$ magnitudes with a Newton-Raphson algorithm, is subdominant. The behavior of the new ISPH scheme has been checked with a number of standard 3D tests that involve large density jumps and contact discontinuities.

A second novelty in the present work has to do with the correct handling of abrupt density jumps. The numerical handling of fluid regions with large density contrasts usually requires us to abandon the Lagrangian formulation of the SPH equations in favor of other schemes. We propose here a self-adaptive scheme, steered by a single parameter of 0 ≤ σ ≤ 1 (see Eqs. (29), (32), and (37)), which selectively chooses the optimal integration scheme. When σ ≃ 0, the Lagrangian-compatible scheme is recovered, while σ ≃ 1 is more adequate to handle large gradients.

The proposed scheme with σ = [0 − 1] gets rid of the tensile instability around contact discontinuities, as it is the closest possible to the compatibility with the Lagrangian formalism because the fraction of particles with a significant σ value is usually low. We can therefore expect to obtain a better depiction of the entropy evolution in shocks than in completely crossed schemes (i.e., with σ = 1). According to Eq. (37), a value of σ ≤ 0.1 is taken in those fluid regions with $\frac{\mathrm{\Delta }\rho }{\overline{\rho }}\le 22\%$, meaning that acoustic waves and weak shock-waves can be handled with, or bordering, the uncrossed option σ ≃ 0. Such moderate density ratios over the kernel domain are also found in self-gravitational structures in equilibrium and in subsonic turbulence experiments.

Table 6 presents a summary of the above-mentioned contributions and contrasts them with our previous works, which sets the current extension into context. Following García-Senz et al. (2012), Cabezón et al. (2012), our chosen method was ISPH, which uses the integral approach to the derivatives, IAD₀, to implement the equations. Further developments were bound to the inclusion of generalized volume elements that improved the partition of the unity.

The results of the hydrodynamic tests unambiguously support the conclusions drawn by the analytical arguments and simple toy models described in Sect. 4. In particular, we find that the proposed volume elements improve both the partition of the unity and the ⟨Δr⟩ condition in all studied cases. Nevertheless, the level of enhancement in the hydrodynamic simulations is lower than that promised by the toy model experiments, especially those concerning the ⟨Δr⟩ constraint. Such a degradation is attributable to the larger particle disorder in the real three-dimensional simulations, but still the results with the improved VEs are valuable. On another note, the novel self-adaptive σ scheme works very well, providing much better results in the hydrostatic and wind-cloud tests than the σ = 0 calculation. The results of the Kelvin–Helmholtz test suggest that the artificial viscosity algorithm plays a central role here, with the VEs and σ choice being subdominant. Finally, the Sedov test clearly shows the best performance of the proposed VEs, which allow us to reproduce the correct density jump across the shock, even in three-dimensional calculations with moderate number of particles. Interestingly, the improved evolution of entropy in shocked regions when the σ = [0 − 1] scheme is used (see Sect. 6.5) reinforces the recent proposal by Rosswog (2020a), who suggested using the entropy as a useful variable for building physically motivated AV triggers in SPH.

As an immediate prospect, we plan to use our improved ISPH code to numerically reproduce the isothermal as well as sub- and supersonic turbulence. Works to implement the ISPH scheme to handle magneto-hydrodynamics effects are also underway.

Acknowledgments

We thank the anonymous referee for insightful comments and suggestions that helped to greatly improve this manuscript. This work has been supported by the MINECO Spanish project PID2020-117252GB-100 (D.G.), by the Swiss Platform for Advanced Scientific Computing (PASC) project SPH-EXA: Optimizing Smooth Particle Hydrodynamics for Exascale Computing (R.C. and D.G.). The authors acknowledge the support of sciCORE (http://scicore.unibas.ch/) scientific computing core facility at University of Basel, where part of these calculations were performed.

References

Appendix A: SPH and ISPH formalisms with generalized VE

According to the Euler–Lagrange formulation of SPH Springel (2010) and references therein) the movement equations are expressed as:

$$\begin{array}{c}\hfill m_a{\ddot{\mathbf{r}}}_a=-{\displaystyle \sum _b{m}_b\frac{P_b}{{\rho }_b^2}\frac{\partial {\rho }_b}{\partial {\mathbf{r}}_a}=\sum _b{P}_b\frac{\partial V_b}{{\partial \mathbf{r}}_a},}\end{array}$$(A.1)

where V_b is the characteristic volume occupied by the particle, and ρ_b = m_b/V_b. The derivative on the RHS embodies the effect of h-gradients. The spatial part of the derivative can be performed in the standard way, ∇W_ab or, better, with the integral approach given by Eq. (31)

We first calculate the value of ∂V_b/∂r_a in Eq. (A.1),

$$\begin{array}{c}\hfill \frac{\partial V_b}{{\partial \mathbf{r}}_a}={\mathrm{\nabla }}_a{V}_b+\frac{\partial V_b}{\partial h_b}\frac{\partial h_b}{\partial {\mathbf{r}}_a}.\end{array}$$(A.2)

The grad-h part of equation above can be estimated differentiating the constraint $h_b^3{V}_b^{-1}=C$ with respect r_a which, after some algebra, gives:

$$\begin{array}{c}\hfill \frac{\partial V_b}{\partial h_b}\frac{\partial h_b}{\partial {\mathbf{r}}_a}[1-\frac{3V_b}{h_b}{\left(\frac{\partial V_b}{\partial h_b}\right)}^{-1}]=-{\mathrm{\nabla }}_a{V}_b,\end{array}$$(A.3)

when combined with expression (A.2) gives:

$$\begin{array}{c}\hfill \frac{\partial V_b}{{\partial \mathbf{r}}_a}={\mathrm{\nabla }}_a{V}_b{[1-\frac{h_b}{3V_b}\frac{\partial V_b}{\partial h_b}]}^{-1},\end{array}$$(A.4)

where ∇_aV_b refers to the spatial gradient whereas the grad-h effects are included in the term in brackets.

To estimate ∇_aV_b, the precise form of the volume elements has to be known. A general form for these elements comes from Hopkins (2013), namely:

$$\begin{array}{c}\hfill V_b=\frac{X_b}{k_b}\mathrm{with}{k}_b={\displaystyle \sum _c{X}_c{W}_{\mathit{bc}}(h_b),}\end{array}$$(A.5)

where X_b is a scalar estimator. Here we have considered two different estimator families, leading to slightly different expressions of the movement and energy equations:

(a) Constant X_b, as for example X_b = m_b and X_b = 1 both reproducing the standard volume elements V_b = m_b/ρ_b. Other choice is $X_b=P_b^k$ where P is the pressure and k ≤ 1. In this case, the estimator is strictly constant only in isobaric systems.

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_a{V}_b\to [b=a]{\mathrm{\nabla }}_a{V}_a& ={\mathrm{\nabla }}_a\left(\frac{X_a}{k_a}\right)=-\frac{X_a}{k_a^2}{\mathrm{\nabla }}_a{k}_a\hfill \\ \hfill & =-\frac{X_a}{k_a^2}{\displaystyle \sum _b{X}_b{\mathrm{\nabla }}_a{W}_{\mathit{ab}}(h_a),}\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_a{V}_b\to [b\ne a]{\mathrm{\nabla }}_a{V}_b& ={\mathrm{\nabla }}_a\left(\frac{X_b}{k_b}\right)=-\frac{X_b}{k_b^2}{\mathrm{\nabla }}_a{k}_b\hfill \\ \hfill & =-\frac{X_b}{k_b^2}{X}_a{\mathrm{\nabla }}_a{W}_{\mathit{ab}}(h_b).\hfill \end{array}$$(A.7)

Combining expressions (A.4), (A.6), and (A.7) with Eq. (A.1), and making use of Eq. (31) to carry out the IA approach of the kernel gradient, the i-component of the acceleration of particle a is finally obtained:

$$\begin{array}{c}\hfill {\ddot{x}}_{i,a}=-\frac{X_a}{m_a}{\displaystyle \sum _b[\frac{X_b{P}_a}{{\mathrm{\Omega }}_a{k}_a^2}{\mathcal{A}}_{i,ab}(h_a)+\frac{X_b{P}_b}{{\mathrm{\Omega }}_b{k}_b^2}{\mathcal{A}}_{i,ab}(h_b)],}\end{array}$$(A.8)

with

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_a=[1-\frac{h_a}{3V_a}\frac{\partial V_a}{\partial h_a}]=[1+\frac{h_a}{3{\rho }_a}\frac{\partial {\rho }_a}{\partial h_a}].\end{array}$$(A.9)

(b) Adaptive $X_b=m_b/{\rho }_b^0$, being ${\rho }_b^0=\sum _c{m}_c{W}_{bc}(h_b)$ the standard SPH density. Because of the normalization condition, this choice leads to a constant k_b = ∑_cX_cW_bc(h_b). Now we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\nabla }}_a{V}_b{\underset{b=a}{\overset{\mathrm{\nabla }}{\to }}}_a{V}_a={\mathrm{\nabla }}_a\left(\frac{X_a}{k_a}\right)=\frac{{\mathrm{\nabla }}_a{X}_a}{k_a}& =-\frac{m_a}{{({\rho }_a^0)}^2{k}_a}{\mathrm{\nabla }}_a{\rho }_a^0\hfill \\ \hfill & =-\frac{m_a}{{({\rho }_a^0)}^2{k}_a}{\displaystyle \sum _b{m}_b{\mathrm{\nabla }}_a{W}_{\mathit{ab}}(h_a),}\hfill \end{array}\end{array}$$(A.10)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\nabla }}_a{V}_b{\underset{b\ne a}{\overset{\mathrm{\nabla }}{\to }}}_a{V}_b={\mathrm{\nabla }}_a\left(\frac{X_b}{k_b}\right)=\frac{{\mathrm{\nabla }}_a{X}_b}{k_b}& =-\frac{m_b}{{({\rho }_b^0)}^2{k}_b}{\mathrm{\nabla }}_a{\rho }_b^0\hfill \\ \hfill & =-\frac{m_b{m}_a}{{({\rho }_b^0)}^2{k}_b}{\mathrm{\nabla }}_a{W}_{\mathit{ab}}(h_b).\hfill \end{array}\end{array}$$(A.11)

Combining expressions (A.4), (A.10), and (A.11) with Eq. (A.1), and making use of Eq. (31) to compute the IA approach of the kernel gradient, the i-component of the acceleration of particle a is obtained,

$$\begin{array}{c}\hfill {\ddot{x}}_{i,a}=-{\displaystyle \sum _b{m}_b[\frac{X_a^2{P}_a}{{\mathrm{\Omega }}_a{m}_a^2{k}_a}{\mathcal{A}}_{i,ab}(h_a)+\frac{X_b^2{P}_b}{{\mathrm{\Omega }}_b{m}_b^2{k}_b}{\mathcal{A}}_{i,ab}(h_b)].}\end{array}$$(A.12)

According to expression (A.9), to estimate Ω_a it is necessary to know the derivative of the density with respect the smoothing-length (∂ρ_a/∂h_a). The result relies in the choice the estimator X_a used to compute the density.

• For constant X_a, as for example (but not necessarily), the standard choice X_a = m_a:

  $$\begin{array}{c}\hfill \frac{\partial {\rho }_a}{\partial h_a}=\frac{m_a}{X_a}{\displaystyle \sum _b}{X}_b\frac{\partial W_{\mathit{ab}}(h_a)}{\partial h_a}.\end{array}$$(A.13)

• The choice $X_a=(m_a/{\rho }_a^0)$ requires a bit more algebra:

  $$\begin{array}{c}\hfill \frac{\partial {\rho }_a}{\partial h_a}=\frac{m_a}{X_a}\left(\frac{\partial k_a}{\partial h_a}\right)-\frac{m_a{k}_a}{X_a^2}\left(\frac{\partial X_a}{\partial h_a}\right),\end{array}$$(A.14)

  with:

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial k_a}{\partial h_a}& ={\displaystyle \sum _b\frac{\partial }{\partial h_a}\left(\frac{m_b}{{\rho }_b^0}\right)W_{\mathit{ab}}(h_a)+\sum _b{X}_b\frac{\partial W_{\mathit{ab}}(h_a)}{\partial h_a}}\hfill \\ \hfill & =-\left(\frac{m_a}{{({\rho }_a^0)}^2}\right)W_{\mathit{aa}}(h_a)\frac{\partial {\rho }_a^0(h_a)}{\partial h_a}+{\displaystyle \sum _b{X}_b\frac{\partial W_{\mathit{ab}}(h_a)}{\partial h_a},}\hfill \end{array}\end{array}$$(A.15)

  and,

  $$\begin{array}{c}\hfill \frac{\partial X_a}{\partial h_a}=\frac{\partial }{\partial h_a}\left(\frac{m_a}{{\rho }_a^0}\right)=-\left(\frac{m_a}{{({\rho }_a^0)}^2}\right)\frac{\partial {\rho }_a^0(h_a)}{\partial h_a}=-\left(\frac{m_a}{{({\rho }_a^0)}^2}\right){\displaystyle \sum _b{m}_b\frac{\partial W_{\mathit{ab}}(h_a)}{\partial h_a}.}\end{array}$$(A.16)

  Setting expressions (A.15) and (A.16) into Eq. (A.14) and via the ensuing manipulation, we have:

  $$\begin{array}{c}\hfill \frac{\partial {\rho }_a}{\partial h_a}=[\frac{{\rho }_a}{{\rho }_a^0}-X_a{W}_{\mathit{aa}}(h_a)]{\displaystyle \sum _b{m}_b\frac{\partial W_{\mathit{ab}}(h_a)}{\partial h_a}+\frac{m_a}{X_a}\sum _b{X}_b\frac{\partial W_{\mathit{ab}}(h_a)}{\partial h_a}.}\end{array}$$(A.17)