Issue 
A&A
Volume 659, March 2022



Article Number  A175  
Number of page(s)  18  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202141877  
Published online  23 March 2022 
Conservative, densitybased smoothed particle hydrodynamics with improved partition of the unity and better estimation of gradients
^{1}
Departament de Física. Universitat Politècnica de Catalunya (UPC), Avinguda Eduard Maristany 16, 08019 Barcelona, Spain
email: domingo.garcia@upc.edu
^{2}
Center for Scientific Computing – sciCORE, Universität Basel, Klingelbergstrasse 61, 4056 Basel, Switzerland
email: ruben.cabezon@unibas.ch
^{3}
Euclid Science Data Center – Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany
email: jescarti@mpe.mpg.de
Received:
26
July
2021
Accepted:
26
December
2021
Context. The accurate evaluation of gradients is a cornerstone of the smoothed particle hydrodynamics (SPH) technique. Using an integral approach to estimating gradients has been proven to substantially enhance its accuracy, retaining the Lagrangian structure of SPH equations and remaining fully conservative. However, in practice, it is difficult to ensure that the Lagrangian formulation is entirely consistent with regard to the exact partition of the unity.
Aims. In this paper, we focus our study on the connection between the choice of the volume elements (VEs) in the SPH summations as well as the accuracy in the gradient estimation within the integral approach scheme (ISPH). We propose a new variant of VEs to improve the partition of the unity that is fully compatible with the Lagrangian formulation of SPH, including gradh corrections.
Aims. Using analytic considerations, simple static toy models in 1D, and a set of full 3D test cases, we show that any improvement in the partition of the unity also leads to a better calculation of gradients when the integral approach is used jointly. Additionally, we propose an easytoimplement modification of the ISPH scheme, which makes it more flexible and better suited to handling sharp density contrasts.
Methods. The ISPH code that is built with the proposed scheme has been validated with a good number of standard tests, some of them involving contact discontinuities. The performance of the code was shown to be excellent in all of these tests, consistently demonstrating that an improvement in the partition of the unity is not detrimental to the optimal conservation of energy, momentum, and entropy that is typical of Lagrangian schemes.
Results. We successfully built a new ISPH scheme on a Lagrangian basis, which is fully conservative, and compatible with selfconsistent gradh terms and an improved partition of the unity. The ensuing code is able to successfully cope with the tensile instability and has been validated with a number of hydrodynamic tests with good results.
Key words: methods: numerical / hydrodynamics
© ESO 2022
1. Introduction
The smoothed particle hydrodynamics (SPH) is a firmly settled numerical technique that is capable of successfully simulating many cuttingedge 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 E0errors, even in the presence of constant pressure fields in nonuniform particle distributions (Agertz et al. 2007; Zhu et al. 2015). One alternative adapts the movingleast 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íaSenz 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_{0} in the seminal paper by GarcíaSenz et al. (2012). It has been shown that the ISPH formulation has the same conservation properties as the standard Lagrangianderived SPH (Springel & Hernquist 2002). We note that a particularly remarkable feature of ISPH is that it effectively reduces the E0error in the derivatives (GarcíaSenz 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 onedimensional 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 reformulated, 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 gradh terms. A second improvement is the use of a novel selfadaptive 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:
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^{1},
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íaSenz et al. 2012):
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,
with ∇f_{a} defined below. The RHS in the integral expression in Eq. (3) becomes:
where
Setting Eqs. (4) and (5) into Eq. (3) allows us to obtain the gradient ∇f_{a},
where
and where the spatial coordinates are represented by x with subindexes 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íaSenz 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^{2}.
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 nonconservative 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,
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:
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:
where 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 GomezGesteira 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 gradh 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 weaklycompressible fluids. Having an almost constant density field for the entire simulation implies that the smoothing length also remains constant and the gradh 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 selfconsistent treatment of the gradh 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^{3} 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:
with k_{a} = ∑_{b}X_{b}W_{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 , 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:
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 timestep, , to estimate the volume elements and density at the current iteration n:
The second is to make use of the density calculated in the standard way in the current timestep, , to estimate X_{a} and then ( calculated with Eq. (12)). That is:
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 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 Lagrangianconsistent scheme (Appendix A), which incorporates the gradh 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_{1}, hereafter) of particle a, as in:
usually improves the requirement:
where E_{2} 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_{1}, E_{2}) will potentially improve the dynamic evolution of the simulated physical system.
4. Estimating the errors E_{1} and E_{2} for different particle distributions
A good control of errors E_{1} and E_{2} is of utmost importance to the SPH technique. This is because the quality of both the interpolated function ⟨f(r)⟩,
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,
Using the Gaussian kernel , the standard SPH derivative of G(x) is expressed as:
thus,
We note that the LHS of Eq. (21) is in fact E_{2}, suggesting that having a good partition of the unity (G ≃ 1, i.e., E_{1} → 0) makes dG/dx ≃ 0 and, as a consequence, the E_{2} error is suppressed. An independent proof of the link between E_{1} and E_{2} 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_{1} and E_{2} can be obtained by studying a handful of static particle distributions and using kernels with compactsupport. These errors are expected to be large close to discontinuities, which in SPH are usually spread over a few times the smoothinglength distance (h). We chose three representative discontinuities that often appear in practical calculations: a Gaussian (model A), an invertedGaussian (model B), and a wall (model C). These are given by the following mathematical expressions:
where the values of the parameters ρ_{0}, Δρ, and δ are specified in Table 1.
Value of the different parameters in profiles A, B, and C that mimic different types of sharp density gradients.
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 sinckernel with exponent n = 5 (Cabezón et al. 2008), which has a shape similar to the M_{6} spline. The different profiles of errors E_{1} and E_{2} at each point, for the estimators X_{1} and X_{2}, are depicted in Fig. 1. As we can see, X_{1} leads to a clear improvement in the kernel normalization condition as p increases. For p ≃ 1 the E_{1} error becomes negligible, hence it is not shown. Interestingly, the estimator X_{2} with p = 1 (black crosses) leads to similar E_{1} and E_{2} errors as X_{1} with p = 0.8 in all profiles. The normalized error E_{2} = ⟨Δ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 (InvertedGaussian), and C (Wall), of Table 1, calculated with X_{1} using p = 0.0 (magenta) and p = 0.8 (blue) in Eq. (13). Crosses (×) in black are for X_{2} (Eq. (15), with p = 1) and points in light blue are the analytic values. Central row: E_{1} (partition of unity) error for the same three profiles following the same color schemes, with the black lines for X_{2} with p = 1. Bottom row: same as central row but for error E_{2} = ⟨Δ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_{1}, E_{2} 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 walllike function f given by profile D in Table 1, Eq. (24)^{4}. Such a derivative, , is sensitive to the choice of the estimator X to compute the VEs. We can thus compare the analytic and the numerical value of and carry out the L_{1} analysis of the results with Eq. (26).
Figure 2 depicts the walllike function f(x) (black solid line) and its analytical gradient of df/dx (lightblue line). Although the value around the maximum of 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_{2} with p = 1 (black dots).
Fig. 2. Walllike function f(x) (black line) and its analytical gradient df/dx (lightblue solid line). The solid lines in magenta and blue were obtained with X_{1} and p = 0.0, 0.9 respectively, while the black points are for X_{2} with p = 1. 
The left panel in Fig. 3 shows the averaged L_{1} value for E_{1} (Eq. (16)) and E_{2} (Eq. (17)):
Fig. 3. L_{1} calculation of the errors in the numerical experiments shown in Fig. 2. Left: averaged value of L_{1} for the partition of the unity (solid lines) and ⟨Δx⟩/h (dotted lines). Blue lines are for the estimator X_{1} and magenta lines for X_{2}. Right: averaged value of L_{1}, obtained with Eq. (26), for the derivative of the wall function calculated with X_{1} (solid blue line) and X_{2} (solid magenta line). The black line is the calculation with the nonconservative 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 X_{1} than for X_{2}. 
calculated in the interval 0.3 ≤ x ≤ 0.7 for the test presented above. The L_{1} 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_{1} or X_{2}. As the figure shows, there is a factor ten reduction of the L_{1} errors for X_{1} and p ≃ 0.9. The errors become negligible when p ≃ 1. Although the option X_{2} shows a lower convergence rate at larger p values, it still provides a significant reduction of L_{1}, decreasing the errors in almost a factor ten when p = 1.
The panel on the right in Fig. 3 shows the averaged L_{1} error of normalized to ,
and how both schemes (ISPH and ncISPH) converge in light of the flavor used for the estimator X (X_{1} solid blue line and X_{2} 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_{1}(p) decreases linearly, achieving the same accuracy as the ncISPH scheme when p > 0.9 for X_{1} or close to it for X_{2}.
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_{2} (Eq. (15)) with p = 1, because, unlike X_{1}, it fits in with the Lagrangian formulation of the SPH equations perfectly and also allows us to eliminate the parameter, p. Additionally, the gradh 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^{5}. 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:
Unfortunately, finding the VEs with such a direct approach generally leads to nonphysical 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 nonphysical. Therefore, we discarded the fully implicit route to finding the VEs in this work.
Fig. 4. Implicit versus explicit estimation of the volume elements, V_{b}, 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:
with V_{a} defined by Eq. (12).
The momentum equation is:
with , 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,
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):
The energy equation is:
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:
where (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, and are the average α and speed of sound between neighboring particles. According to Monaghan (1997), the AV contribution to the energy equation should include a heatconduction term which smooths the pressure in wallshock conditions. The precise form of such heatconduction term is:
with , which in the SPHYNX code is implemented as:
The signal velocity used in our tests is (Price 2008):
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 densitypeak 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 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)^{6} 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^{7}. In spite of not being totally Lagrangiancompatible and with the gradh 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 selfadaptive 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.
with and 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 wellknown 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 twophase 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, highdensity spherical cloud (the windcloud test), as well as to simulate the evolution of a pointlike explosion (the Sedov test). The simulations were carried out using the ISPH code SPHYNX and the outcomes have been compared to wellknown 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^{8} of X_{a} = m_{a} and the improved VEs with .
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^{−6} from the initial value. The relative deviation, ΔE/E_{0} in the total energy was less than 10^{−5}, except in the isobaric twofluid test (Sect. 6.1), where it was ≃10^{−4} at late times.
6.1. Isobaric twofluid test
The simulation of the hydrostatic evolution of a twophase 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:
The system is isobaric with P = 2.5 and we use N = 110^{3} equalmass particles spread in two nested bodycentred 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 nonphysically 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 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_{5}, H_{6}), calculated with σ = 0, the system completely looses its shape in half of the soundcrossing time, τ_{sc} ≃ 0.9 s. Nevertheless, increasing the conductivity in the AV (model H_{6}) slightly delays the deformation; although, in the end, it is not able to prevent it. Model H_{7}, 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 densitybased 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_{8}, 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_{7}. As shown in Fig. 6 (upper panel), the algorithm (Eq. (37)) that selfadapts σ 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 H_{5}, H_{6}, H_{7}, H_{8} (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 soundcrossing time being τ_{sc} = 0.9. 
Fig. 6. Isobaric twofluid numerical experiment. Upper panel: slice around Z = 0 showing the averaged σparameter (Eq. (37)), in model H_{8} at t = 0.55. Appreciable values of σ are only attained at the fluid interphase. Lower panel: time evolution of averaged L_{1}(E_{1}) (solid lines), L_{1}(E_{2}) (dashed lines) for models H_{3}, H_{7}, and H_{8} in Table 2. 
L_{1} values for errors E_{1} and E_{2} at t = 1.5 in the hydrostatic square test.
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_{1} values (Eq. (25)) in a rectangular shell with width 0.3 around the contact discontinuity.
The column L_{1}(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:
where x, y, and z are the Cartesian coordinates of the particles while x_{0}, y_{0}, and z_{0} coordinates represent their initial position. Subindex b runs from 1 up to the N particles contained in the rectangular shell.
The full Lagrangian models H_{1}, H_{2}, H_{5}, and H_{6}, calculated with σ = 0, give the worst results, with L_{1}(d)≥3 × 10^{−2} at t = 1.5. The minimum value of L_{1}(d) corresponds to model H_{7} (L_{1}(d) = 1.5 × 10^{−2}), calculated with σ = 1 and improved VEs. However, model H_{8} calculated with the selfadaptive σ[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_{1} of the partition of the unity (solid lines) and ⟨Δr⟩/h (dashed lines) for models H_{3}, H_{7}, and H_{8} in Table 2. Comparing H_{3} to H_{7} highlights the relevance of the VEs, while comparing H_{7} to H_{8} shows the impact of a variable σ. Is it clear that the partition of the unity is enhanced with the choice X = m/ρ^{0}, nevertheless the behavior of ⟨Δr⟩/h is less defined. The error L_{1}(E_{2}) remains lower for models H_{7} and H_{8} than model H_{3} until t ≃ 0.2. Afterwards the error grows to eventually reach a similar value to that of model H_{3}, 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/ρ^{0} 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íaSenz 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 threedimensional layer of size [1 × 1 × 0.0625] with 4.2 × 10^{6} equalmass 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 highdensity one. Each region was generated from a random particle distribution relaxed to a glasslike 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 ) 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.
Simulated models in the Kelvin–Helmholz test.
In Fig. 7, we represent the particle distributions for each simulated model (each snapshot corresponds to one model, from KH_{1} to KH_{4}) 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 KH_{1} to KH_{4} at t = 2 (t_{KH} = 1.06). The density [1:2] is colorcoded. 
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 v_{y} 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_{1} value for errors E_{1} and E_{2}. To average the L_{1} values we restricted to those particles that had σ ≥ 0.025. This condition returns all those particles that are found in the interphase (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_{1} for all simulated tests at t = 2, where we can see that the choice X = m/ρ^{0} reduces both, the L_{1}(E_{1}) and the L_{1}(E_{2}) errors.
Fig. 9. L_{1} evolution of errors E_{1} and E_{2} for the KH models. Solid lines correspond to X = m, while dashed lines are those of X = m/ρ^{0}. 
From these results, we can conclude that despite decreasing the errors E_{1} and E_{2}, 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_{5} and KH_{6}, 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_{5} and KH_{6} 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 lowdensity regions by a factor of 4, simply by generating a new relaxed model with four times fewer particles, along with repeated simulations KH_{2} and KH_{4}, now with this new density jump by a factor of 8 (simulations KH_{7} and KH_{8}). 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 nonlinear 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 subdominant. 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. Windcloud collision
The windcloud 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 twophase medium with a large density contrast. In this test, a spherical cloud of cold gas, initially at rest, is swept by a lowdensity 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 windcloud 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^{6} and N_{c} = 1.23 × 10^{5} particles for the wind and the cloud, respectively. The initial particle distribution for both was that of a stable glasslike configuration. The box has a size {0.25 × 0.25 × 1} with periodic boundary conditions. The highdensity 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_{1} values of errors E_{1} and E_{2} 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_{1} 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_{1} error generated by the particles in the interphase. 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. L_{1} values for E_{1} (solid lines) and E_{2} (dashed lines) in the windcloud collision tests. Using the improved VEs reduces both L_{1} values for the whole simulation. 
Fig. 12. Particle distribution of the windcloud test in a thin slice around Y = 0 for model W_{6} at t/t_{KH} = 0, 1, 2, and 3. The parameter σ is colorcoded. 
Simulated models in the windcloud test.
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 shockwave is born, which moves through the dense cloud and compresses its matter. The density jumps from the original value of ρ_{0} = 10 and tends toward the limiting value ρ = 4ρ_{0}, 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/ρ^{0}, give in general larger values than models calculated with X = m.
Fig. 13. Particle distribution of the windcloud test in a thin slice around Y = 0 for models W_{3} and W_{6} at t_{KH} = 1.5. Density is colorcoded. 
In order to characterize the destruction of the highdensity 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/ρ^{0} 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/t_{KH} for the windcloud 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 threedimensional square box of side L = 1 and total number of particles N = 100^{3} arranged in a glasslike 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_{1} errors in the shocked region^{9} 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_{4} and S_{8}. In both cases, the selfadaptive algorithm to estimate σ was active. The best match with the analytic profile (black line) was obtained with , 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_{4}) and improved VEs (S_{8}), respectively. Figure 17 shows the distribution of the averaged σ in a thin slice along the system for model S_{8}. As it can be seen, the value of σ selfadapts so that it approaches to 1 across the shock front, goes down to ≃0.5 in the postshock region and vanishes in the central region.
Fig. 16. Sedov explosion for models S_{4} and S_{8} in Table 5. Top row: density profiles at t = 0.09 calculated with X_{a} = m_{a} (model S_{4}, standard VEs) and (model S_{8}, 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 S_{8} in Table 5 at t = 0.09. 
Fig. 18. Sedov explosion. Evolution of the averaged L_{1} error, (Eq. (25)), in the shocked region corresponding to models S_{4}, (X_{a} = m_{a}), and S_{8}, (), in Table 5. We show both the partition of unity and the normalized ⟨Δr⟩/h condition, for the estimator choices: X_{a} = m_{a} (purple and orange lines) and (blue and magenta lines). Also shown is the evolution of the maximum density for the two volume estimators (dashed lines, black is for S_{4} and lightblue for S_{8}) during the explosion. 
Summary of all calculated Sedov tests stating the values of σ, heat transport constant in the AV, VE choice, density peak, and L_{1} errors at t = 0.09.
The profile of the normalization function ∑_{b}V_{b}W_{ab}(h_{a}) at t = 0.09 is also shown in Fig. 16, bottom row. As it can be seen, the choice 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_{1}(E_{1}, E_{2}) 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 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_{1} to ρ_{max}(t = 0.09) = 3.16 in model S_{2}. 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 Lagrangiancompatible. 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:
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:
with u = c_{v}T. Thus, integrating on both sides:
where ρ_{0}, u_{0} 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:
which, integrating on both sides, leads to:
where 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):
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 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 Lagrangiancompatible estimator, such as , 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 Lagrangiancompatibility and the suppression of the tensile instability.
Fig. 19. Radial profile of the entropy variation by pressure forces, Δs_{P}/c_{v}, for weak (upper row, at t = 0.18) and strong (lower row, at t = 0.08) pointlike explosions and different σ choices. The first column is for estimator X_{a} = m_{a} and the second for . 
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, Δs_{tot}/c_{v}, 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 postshock 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 ∑_{b}V_{b}W_{ab}(h_{a}) = 1 automatically leads to an enhancement of the condition ⟨Δr⟩=∑_{b}V_{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 gradh terms in the momentum and energy equations. To do that the volume elements, V_{a}, have been redefined so that V_{a} = X_{a}/k_{a} with ; and . That particular choice of V_{a} results in the new set of ISPH equations described in Appendix A. These ISPH equations are not only Lagrangiancompatible, 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 magnitudes with a NewtonRaphson 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 selfadaptive 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 Lagrangiancompatible 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 , meaning that acoustic waves and weak shockwaves can be handled with, or bordering, the uncrossed option σ ≃ 0. Such moderate density ratios over the kernel domain are also found in selfgravitational structures in equilibrium and in subsonic turbulence experiments.
Table 6 presents a summary of the abovementioned contributions and contrasts them with our previous works, which sets the current extension into context. Following GarcíaSenz et al. (2012), Cabezón et al. (2012), our chosen method was ISPH, which uses the integral approach to the derivatives, IAD_{0}, to implement the equations. Further developments were bound to the inclusion of generalized volume elements that improved the partition of the unity.
Summary of the features of the approach presented in this work within the context of the works on which it is based.
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 threedimensional simulations, but still the results with the improved VEs are valuable. On another note, the novel selfadaptive σ scheme works very well, providing much better results in the hydrostatic and windcloud 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 threedimensional 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 magnetohydrodynamics effects are also underway.
Admittedly, this is not the most used implementation of the derivative in SPH. The most commonly used version is a variation involving differences of the derived function. Nevertheless, it is formally the same as Eq. (2).
The widely used Intel MKL PARDISO library function is based on a legacy version of this project: https://www.pardisoproject.org/
The possibility of working with (with some preference for σ = 0) was already suggested by Monaghan (1992) in his popular review, but σ was thought to be a constant real number.
Specifically, it can be shown that the scheme in Paper I involve terms such as , with p = 0.7 being the exponent of the explicit estimator (Eq. (14)) considered in that work.
Actually, any constant magnitude is a suitable choice as X_{a}, but taking the mass of the particle allows for the finetuning of the density (if needed) by slightly modifying the mass of the particles. In nearly isobaric systems, with shallow pressure gradients, the choice could also be appropriate.
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 PID2020117252GB100 (D.G.), by the Swiss Platform for Advanced Scientific Computing (PASC) project SPHEXA: 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
 Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963 [NASA ADS] [CrossRef] [Google Scholar]
 Bonet, J. 1999, Comput. Methods Appl. Mech. Eng., 180, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Brookshaw, L. 1985, PASA, 6, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Cabezón, R. M., GarcíaSenz, D., & Relaño, A. 2008, J. Comput. Phys., 227, 8523 [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]
 Cabezón, R. M., GarcíaSenz, D., & Figueira, J. 2017, A&A, 606, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Colagrossi, A., & Landrini, M. 2003, J. Comput. Phys., 191, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669 [NASA ADS] [CrossRef] [Google Scholar]
 Dilts, G. A. 1999, Int. J. Numer. Methods Eng., 44, 1115 [NASA ADS] [CrossRef] [Google Scholar]
 Escartín, J. 2016, Ph.D. Thesis, Polytechnical University of Catalonia, http://hdl.handle.net/10803/384002 [Google Scholar]
 Frontiere, N., Raskin, C. D., & Owen, J. M. 2017, J. Comput. Phys., 332, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Fulk, D. A., & Quinn, D. W. 1996, J. Comput. Phys., 126, 165 [NASA ADS] [CrossRef] [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 [NASA ADS] [CrossRef] [Google Scholar]
 GomezGesteira, M., Rogers, B., Dalrymple, R., & Crespo, A. J. 2010, J. Hydraulic Res., 48, 6 [CrossRef] [Google Scholar]
 Hopkins, P. F. 2013, MNRAS, 428, 2840 [NASA ADS] [CrossRef] [Google Scholar]
 Knapp, C. E. 2000, Ph.D. Thesis, Los Alamos National Laboratory, https://www.osti.gov/biblio/754046 [Google Scholar]
 Liu, M., & Liu, G. 2006, Appl. Numer. Math., 56, 19 [CrossRef] [Google Scholar]
 Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 McNally, C. P., Lyra, W., & Passy, J.C. 2012, ApJS, 201, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Monaghan, J. J. 1992, ARA&A, 30, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Monaghan, J. J. 1997, J. Comput. Phys., 136, 298 [NASA ADS] [CrossRef] [Google Scholar]
 Monaghan, J. J. 2005, Rep. Prog. Phys., 68, 1703 [NASA ADS] [CrossRef] [Google Scholar]
 Monaghan, J. J. 2012, Ann. Rev. Fluid Mech., 44, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Oger, G., Doring, M., Alessandrini, B., & Ferrant, P. 2007, J. Comput. Phys., 225, 1472 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. 2004, Ph.D. Thesis, University of Cambridge, UK, [arXiv:astroph/0507472] [Google Scholar]
 Price, D. J. 2008, J. Comput. Phys., 227, 10040 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031 [Google Scholar]
 Read, J. I., Hayfield, T., & Agertz, O. 2010, MNRAS, 405, 1513 [NASA ADS] [Google Scholar]
 Ritchie, B. W., & Thomas, P. A. 2001, MNRAS, 323, 743 [NASA ADS] [CrossRef] [Google Scholar]
 Rosswog, S. 2015a, MNRAS, 448, 3628 [NASA ADS] [CrossRef] [Google Scholar]
 Rosswog, S. 2015b, Liv. Rev. Comput. Astrophys., 1, 1 [CrossRef] [Google Scholar]
 Rosswog, S. 2020a, ApJ, 898, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Rosswog, S. 2020b, MNRAS, 498, 4230 [NASA ADS] [CrossRef] [Google Scholar]
 Saitoh, T. R., & Makino, J. 2013, ApJ, 768, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2010, ARA&A, 48, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649 [NASA ADS] [CrossRef] [Google Scholar]
 Tricco, T. S. 2019, MNRAS, 488, 5210 [NASA ADS] [CrossRef] [Google Scholar]
 Valdarnini, R. 2016, ApJ, 831, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Wadsley, J. W., Keller, B. W., & Quinn, T. R. 2017, MNRAS, 471, 2357 [NASA ADS] [CrossRef] [Google Scholar]
 Wissing, R., & Shen, S. 2020, A&A, 638, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhu, Q., Hernquist, L., & Li, Y. 2015, ApJ, 800, 6 [NASA ADS] [CrossRef] [Google Scholar]
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:
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 hgradients. 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),
The gradh part of equation above can be estimated differentiating the constraint with respect r_{a} which, after some algebra, gives:
when combined with expression (A.2) gives:
where ∇_{a}V_{b} refers to the spatial gradient whereas the gradh effects are included in the term in brackets.
To estimate ∇_{a}V_{b}, the precise form of the volume elements has to be known. A general form for these elements comes from Hopkins (2013), namely:
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 where P is the pressure and k ≤ 1. In this case, the estimator is strictly constant only in isobaric systems.
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 icomponent of the acceleration of particle a is finally obtained:
with
(b) Adaptive , being the standard SPH density. Because of the normalization condition, this choice leads to a constant k_{b} = ∑_{c}X_{c}W_{bc}(h_{b}). Now we have
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 icomponent of the acceleration of particle a is obtained,
According to expression (A.9), to estimate Ω_{a} it is necessary to know the derivative of the density with respect the smoothinglength (∂ρ_{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}:

The choice requires a bit more algebra:
with:
and,
Setting expressions (A.15) and (A.16) into Eq. (A.14) and via the ensuing manipulation, we have:
All Tables
Value of the different parameters in profiles A, B, and C that mimic different types of sharp density gradients.
L_{1} values for errors E_{1} and E_{2} at t = 1.5 in the hydrostatic square test.
Summary of all calculated Sedov tests stating the values of σ, heat transport constant in the AV, VE choice, density peak, and L_{1} errors at t = 0.09.
Summary of the features of the approach presented in this work within the context of the works on which it is based.
All Figures
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 (InvertedGaussian), and C (Wall), of Table 1, calculated with X_{1} using p = 0.0 (magenta) and p = 0.8 (blue) in Eq. (13). Crosses (×) in black are for X_{2} (Eq. (15), with p = 1) and points in light blue are the analytic values. Central row: E_{1} (partition of unity) error for the same three profiles following the same color schemes, with the black lines for X_{2} with p = 1. Bottom row: same as central row but for error E_{2} = ⟨Δx⟩/h. 

In the text 
Fig. 2. Walllike function f(x) (black line) and its analytical gradient df/dx (lightblue solid line). The solid lines in magenta and blue were obtained with X_{1} and p = 0.0, 0.9 respectively, while the black points are for X_{2} with p = 1. 

In the text 
Fig. 3. L_{1} calculation of the errors in the numerical experiments shown in Fig. 2. Left: averaged value of L_{1} for the partition of the unity (solid lines) and ⟨Δx⟩/h (dotted lines). Blue lines are for the estimator X_{1} and magenta lines for X_{2}. Right: averaged value of L_{1}, obtained with Eq. (26), for the derivative of the wall function calculated with X_{1} (solid blue line) and X_{2} (solid magenta line). The black line is the calculation with the nonconservative 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 X_{1} than for X_{2}. 

In the text 
Fig. 4. Implicit versus explicit estimation of the volume elements, V_{b}, in the numerical experiments with the Gaussian curve (Eq. (22)). 

In the text 
Fig. 5. Slices around Z = 0 showing the density colormap of models H_{5}, H_{6}, H_{7}, H_{8} (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 soundcrossing time being τ_{sc} = 0.9. 

In the text 
Fig. 6. Isobaric twofluid numerical experiment. Upper panel: slice around Z = 0 showing the averaged σparameter (Eq. (37)), in model H_{8} at t = 0.55. Appreciable values of σ are only attained at the fluid interphase. Lower panel: time evolution of averaged L_{1}(E_{1}) (solid lines), L_{1}(E_{2}) (dashed lines) for models H_{3}, H_{7}, and H_{8} in Table 2. 

In the text 
Fig. 7. Particles distribution in a thin slice around Z = 0 for models KH_{1} to KH_{4} at t = 2 (t_{KH} = 1.06). The density [1:2] is colorcoded. 

In the text 
Fig. 8. Amplitude growth of the v_{y} 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. 

In the text 
Fig. 9. L_{1} evolution of errors E_{1} and E_{2} for the KH models. Solid lines correspond to X = m, while dashed lines are those of X = m/ρ^{0}. 

In the text 
Fig. 10. Color map of σ (Eq. (37)), in a thin slice around Z = 0 for model KH_{4} at t = 2. 

In the text 
Fig. 11. L_{1} values for E_{1} (solid lines) and E_{2} (dashed lines) in the windcloud collision tests. Using the improved VEs reduces both L_{1} values for the whole simulation. 

In the text 
Fig. 12. Particle distribution of the windcloud test in a thin slice around Y = 0 for model W_{6} at t/t_{KH} = 0, 1, 2, and 3. The parameter σ is colorcoded. 

In the text 
Fig. 13. Particle distribution of the windcloud test in a thin slice around Y = 0 for models W_{3} and W_{6} at t_{KH} = 1.5. Density is colorcoded. 

In the text 
Fig. 14. Evolution of the maximum density in the cloud at t/t_{KH} ≤ 1.6 for all models in Table 4. 

In the text 
Fig. 15. Percent of surviving cloud in function of t/t_{KH} for the windcloud collision tests. 

In the text 
Fig. 16. Sedov explosion for models S_{4} and S_{8} in Table 5. Top row: density profiles at t = 0.09 calculated with X_{a} = m_{a} (model S_{4}, standard VEs) and (model S_{8}, 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. 

In the text 
Fig. 17. Color map slice around Z = 0 showing the value of the averaged σparameter defined in Eq. (37), for model S_{8} in Table 5 at t = 0.09. 

In the text 
Fig. 18. Sedov explosion. Evolution of the averaged L_{1} error, (Eq. (25)), in the shocked region corresponding to models S_{4}, (X_{a} = m_{a}), and S_{8}, (), in Table 5. We show both the partition of unity and the normalized ⟨Δr⟩/h condition, for the estimator choices: X_{a} = m_{a} (purple and orange lines) and (blue and magenta lines). Also shown is the evolution of the maximum density for the two volume estimators (dashed lines, black is for S_{4} and lightblue for S_{8}) during the explosion. 

In the text 
Fig. 19. Radial profile of the entropy variation by pressure forces, Δs_{P}/c_{v}, for weak (upper row, at t = 0.18) and strong (lower row, at t = 0.08) pointlike explosions and different σ choices. The first column is for estimator X_{a} = m_{a} and the second for . 

In the text 
Fig. 20. Radial distribution function of total entropy, Δs_{tot}/c_{v}, 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. 

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.