Free Access
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/0004-6361/202141877
Published online 23 March 2022

© 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 IAD0 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:

(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/ρ. Specifically1,

(2)

where a and b refer to neighbouring particles with masses ma and mb, 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):

(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,

(4)

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

(5)

where

(6)

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

(7)

where

(8)

and where the spatial coordinates are represented by x with sub-indexes i or j. Henceforth, we note Wab(ha)≡W(|rb − ra|,ha) 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 ⟨Δra → 0. We refer to ISPH as the conservative scheme which simply neglects the term f(ra)⟨Δra in Eq. (5). This is justified because ⟨Δr⟩ is, in fact, the indefinite integral of an odd function, which is zero2.

On the other hand, the complete integral approach, which takes into account the f(ra)⟨Δra 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 ⟨Δra ≃ 0. Having both, a perfect partition of unity and ⟨Δra = 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 Va = ma/ρa, which leads to the widely used density equation,

(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:

(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:

(11)

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 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, Va, may be of interest to address specific problems. The code SPHYNX3 from Cabezón et al. (2017), makes use of the concept of generalized volume elements (Hopkins 2013; Saitoh & Makino 2013). First, a scalar estimator Xa is defined so that the particle volume is:

(12)

with ka = ∑bXbWab(ha). The density of the particle is then calculated as ρa = ma/Va. 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:

(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, , to estimate the volume elements and density at the current iteration n:

(14)

The second is to make use of the density calculated in the standard way in the current time-step, , to estimate Xa and then ( calculated with Eq. (12)). That is:

(15)

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

In the simple, static, toy models discussed below, X1, 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, X1, 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 X1, 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 Xa which, in the end, was the one chosen in this work. Therefore, Xa consistently refers to X2, 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 X1, a leads to a better partition of the unity and interpolations than X2, 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 X2, 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 (E1, hereafter) of particle a, as in:

(16)

usually improves the requirement:

(17)

where E2 is the normalized error module |Δr⟩|a/ha of particle a. The error reduction is the consequence of using the estimators X1, 2 to evaluate the volume elements. Any reduction in both errors (E1, E2) will potentially improve the dynamic evolution of the simulated physical system.

4. Estimating the errors E1 and E2 for different particle distributions

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

(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,

(19)

Using the Gaussian kernel , the standard SPH derivative of G(x) is expressed as:

(20)

thus,

(21)

We note that the LHS of Eq. (21) is in fact E2, suggesting that having a good partition of the unity (G ≃ 1, i.e., E1 → 0) makes dG/dx ≃ 0 and, as a consequence, the E2 error is suppressed. An independent proof of the link between E1 and E2 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 E1 and E2 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:

(22)

(23)

(24)

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

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 sinc-kernel with exponent n = 5 (Cabezón et al. 2008), which has a shape similar to the M6 spline. The different profiles of errors E1 and E2 at each point, for the estimators X1 and X2, are depicted in Fig. 1. As we can see, X1 leads to a clear improvement in the kernel normalization condition as p increases. For p ≃ 1 the E1 error becomes negligible, hence it is not shown. Interestingly, the estimator X2 with p = 1 (black crosses) leads to similar E1 and E2 errors as X1 with p = 0.8 in all profiles. The normalized error E2 = |⟨Δ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.

thumbnail 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 E1, E2 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)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 L1 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 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 X2 with p = 1 (black dots).

thumbnail 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 L1 value for E1 (Eq. (16)) and E2 (Eq. (17)):

(25)

thumbnail 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 L1 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, X1 or X2. As the figure shows, there is a factor ten reduction of the L1 errors for X1 and p ≃ 0.9. The errors become negligible when p ≃ 1. Although the option X2 shows a lower convergence rate at larger p values, it still provides a significant reduction of L1, decreasing the errors in almost a factor ten when p = 1.

The panel on the right in Fig. 3 shows the averaged L1 error of normalized to ,

(26)

and how both schemes (ISPH and ncISPH) converge in light of the flavor used for the estimator X (X1 solid blue line and X2 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 L1(p) decreases linearly, achieving the same accuracy as the ncISPH scheme when p > 0.9 for X1 or close to it for X2.

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 X2 (Eq. (15)) with p = 1, because, unlike X1, 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 library5. These types of libraries could also be used to implicitly find the volume elements Vb by solving large linear systems with N × N equations and unknowns:

(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.

thumbnail 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., Xa = ma) 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:

(28)

with Va defined by Eq. (12).

The momentum equation is:

(29)

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,

(30)

with cij, 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):

(31)

The energy equation is:

(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:

(33)

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 heat-conduction term which smooths the pressure in wall-shock conditions. The precise form of such heat-conduction term is:

(34)

with , which in the SPHYNX code is implemented as:

(35)

The signal velocity used in our tests is (Price 2008):

(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 to [Pa/(ρ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. ub/ρa in Ritchie & Thomas 2001 or Pa/(ρaρb) in Read et al. 2010). In Cabezón et al. (2017), the choice of X1, a (Eq. (14)), combined to the SPH equations defined in Rosswog (2015b), naturally led to an almost “crossed” scheme7. 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.

(37)

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 Atmin = 0.1, Atmax = 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 choice8 of Xa = ma 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|/E0 in the total energy was less than 10−5, except in the isobaric two-fluid test (Sect. 6.1), where it was ≃10−4 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:

(38)

The system is isobaric with P = 2.5 and we use N = 1103 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 nb = 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 are shown in Fig. 5, and calculations with Xa = ma 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 H5, H6), 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 H6) slightly delays the deformation; although, in the end, it is not able to prevent it. Model H7, 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 X1, a given by Eq. (14) with p = 0.7 (see footnote 7). It is worth noting that model H8, 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 H7. 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.

thumbnail 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.

thumbnail 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.

Table 2.

L1 values for errors E1 and E2 at t = 1.5 in the hydrostatic square test.

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

The column L1(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:

(39)

where x, y, and z are the Cartesian coordinates of the particles while x0, y0, and z0 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 H1, H2, H5, and H6, calculated with σ = 0, give the worst results, with L1(d)≥3 × 10−2 at t = 1.5. The minimum value of L1(d) corresponds to model H7 (L1(d) = 1.5 × 10−2), calculated with σ = 1 and improved VEs. However, model H8 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 L1 of the partition of the unity (solid lines) and ⟨Δr⟩/h (dashed lines) for models H3, H7, and H8 in Table 2. Comparing H3 to H7 highlights the relevance of the VEs, while comparing H7 to H8 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 L1(E2) remains lower for models H7 and H8 than model H3 until t ≃ 0.2. Afterwards the error grows to eventually reach a similar value to that of model H3, calculated with Xa = ma. 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í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 × 106 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 tKH, as defined in Agertz et al. (2007), which results in tKH = 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 Xa = ma 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.

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 KH1 to KH4) 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.

thumbnail 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 vy 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.

thumbnail 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 L1 value for errors E1 and E2. To average the L1 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 L1 for all simulated tests at t = 2, where we can see that the choice X = m/ρ0 reduces both, the L1(E1) and the L1(E2) errors.

thumbnail 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.

thumbnail 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 E1 and E2, 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, KH5 and KH6, 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 KH5 and KH6 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 vy 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 KH2 and KH4, now with this new density jump by a factor of 8 (simulations KH7 and KH8). We present the evolution of the amplitude growth of the vy 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 Nw = 11.23 × 106 and Nc = 1.23 × 105 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 uw = 3/2 and uc = 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 L1 values of errors E1 and E2 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 L1 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 L1 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.

thumbnail 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.

thumbnail 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.

Table 4.

Simulated models in the wind-cloud test.

The density color map at t/tKH = 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 ρ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.

thumbnail 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.

thumbnail 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.9uw) in function of t/tKH 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).

thumbnail 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 = 1003 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 E1, 2 (Eqs. (16) and (17)) with respect the VEs estimator choice and σ, we calculated the L1 errors in the shocked region9 with Eq. (25).

The outcome of the simulations is summarized in Table 5 and Figs. 1618. 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 S4 and S8. In both cases, the self-adaptive 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 (S4) and improved VEs (S8), respectively. Figure 17 shows the distribution of the averaged σ in a thin slice along the system for model S8. 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.

thumbnail 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 (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.

thumbnail 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.

thumbnail Fig. 18.

Sedov explosion. Evolution of the averaged L1 error, (Eq. (25)), in the shocked region corresponding to models S4,  (Xa = ma), and S8,  (), 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 (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.

Table 5.

Summary of all calculated Sedov tests stating the values of σ, heat transport constant in the AV, VE choice, density peak, and L1 errors at t = 0.09.

The profile of the normalization function ∑bVbWab(ha) 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 Xa = ma. Such assessment is quantitatively confirmed by the temporal evolution of the error estimators L1(E1, E2) 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 S1 to ρmax(t = 0.09) = 3.16 in model S2. 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 ΔsP and that from the artificial viscosity ΔsAV. Ideally, ΔsP ≃ 0 because the dissipation is basically driven by the AV.

The second law of thermodynamics states:

(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, cv, we have:

(41)

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

(42)

where ρ0, u0 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/cv due to the AV. Subtracting (Δs)AV/cv 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:

(43)

which, integrating on both sides, leads to:

(44)

where is directly obtained from the SPH code. The variation of the entropy of the gas, due to pressure forces (Δs)P/cv, and excluding the AV contribution, is determined subtracting Eq. (44) from Eq. (42):

(45)

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

In Fig. 19, we show the radial profile of |ΔsP|/cv 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 Lagrangian-compatible 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 Lagrangian-compatibility and the suppression of the tensile instability.

thumbnail 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 .

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.

thumbnail 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 ∑bVbWab(ha) = 1 automatically leads to an enhancement of the condition ⟨Δr⟩=∑bVb(rb − ra)Wab(ha) = 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, Va, have been re-defined so that Va = Xa/ka with ; and . That particular choice of Va 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 ha and 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 , 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, IAD0, to implement the equations. Further developments were bound to the inclusion of generalized volume elements that improved the partition of the unity.

Table 6.

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 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.


1

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).

2

Interestingly, this argument has been used for the whole history of SPH to justify that it is a second-order method, which is true in the continuum limit, but not always numerically ensured.

4

This test would mimic, for example, a thermal wave passing through a star.

5

The widely used Intel MKL PARDISO library function is based on a legacy version of this project: https://www.pardiso-project.org/

6

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.

7

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.

8

Actually, any constant magnitude is a suitable choice as Xa, but taking the mass of the particle allows for the fine-tuning 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.

9

Defined as the volume with specific internal energy u(t)≥0.1.

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

  1. Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bonet, J. 1999, Comput. Methods Appl. Mech. Eng., 180, 97 [NASA ADS] [CrossRef] [Google Scholar]
  3. Brookshaw, L. 1985, PASA, 6, 207 [NASA ADS] [CrossRef] [Google Scholar]
  4. Cabezón, R. M., García-Senz, D., & Relaño, A. 2008, J. Comput. Phys., 227, 8523 [CrossRef] [Google Scholar]
  5. Cabezón, R. M., García-Senz, D., & Escartín, J. A. 2012, A&A, 545, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Cabezón, R. M., García-Senz, D., & Figueira, J. 2017, A&A, 606, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Colagrossi, A., & Landrini, M. 2003, J. Comput. Phys., 191, 448 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669 [NASA ADS] [CrossRef] [Google Scholar]
  9. Dilts, G. A. 1999, Int. J. Numer. Methods Eng., 44, 1115 [NASA ADS] [CrossRef] [Google Scholar]
  10. Escartín, J. 2016, Ph.D. Thesis, Polytechnical University of Catalonia, http://hdl.handle.net/10803/384002 [Google Scholar]
  11. Frontiere, N., Raskin, C. D., & Owen, J. M. 2017, J. Comput. Phys., 332, 160 [NASA ADS] [CrossRef] [Google Scholar]
  12. Fulk, D. A., & Quinn, D. W. 1996, J. Comput. Phys., 126, 165 [NASA ADS] [CrossRef] [Google Scholar]
  13. García-Senz, D., Cabezón, R. M., & Escartín, J. A. 2012, A&A, 538, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375 [NASA ADS] [CrossRef] [Google Scholar]
  15. Gomez-Gesteira, M., Rogers, B., Dalrymple, R., & Crespo, A. J. 2010, J. Hydraulic Res., 48, 6 [CrossRef] [Google Scholar]
  16. Hopkins, P. F. 2013, MNRAS, 428, 2840 [NASA ADS] [CrossRef] [Google Scholar]
  17. Knapp, C. E. 2000, Ph.D. Thesis, Los Alamos National Laboratory, https://www.osti.gov/biblio/754046 [Google Scholar]
  18. Liu, M., & Liu, G. 2006, Appl. Numer. Math., 56, 19 [CrossRef] [Google Scholar]
  19. Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] [Google Scholar]
  20. McNally, C. P., Lyra, W., & Passy, J.-C. 2012, ApJS, 201, 18 [NASA ADS] [CrossRef] [Google Scholar]
  21. Monaghan, J. J. 1992, ARA&A, 30, 543 [NASA ADS] [CrossRef] [Google Scholar]
  22. Monaghan, J. J. 1997, J. Comput. Phys., 136, 298 [NASA ADS] [CrossRef] [Google Scholar]
  23. Monaghan, J. J. 2005, Rep. Prog. Phys., 68, 1703 [NASA ADS] [CrossRef] [Google Scholar]
  24. Monaghan, J. J. 2012, Ann. Rev. Fluid Mech., 44, 323 [NASA ADS] [CrossRef] [Google Scholar]
  25. Oger, G., Doring, M., Alessandrini, B., & Ferrant, P. 2007, J. Comput. Phys., 225, 1472 [NASA ADS] [CrossRef] [Google Scholar]
  26. Price, D. 2004, Ph.D. Thesis, University of Cambridge, UK, [arXiv:astroph/0507472] [Google Scholar]
  27. Price, D. J. 2008, J. Comput. Phys., 227, 10040 [NASA ADS] [CrossRef] [Google Scholar]
  28. Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031 [Google Scholar]
  29. Read, J. I., Hayfield, T., & Agertz, O. 2010, MNRAS, 405, 1513 [NASA ADS] [Google Scholar]
  30. Ritchie, B. W., & Thomas, P. A. 2001, MNRAS, 323, 743 [NASA ADS] [CrossRef] [Google Scholar]
  31. Rosswog, S. 2015a, MNRAS, 448, 3628 [NASA ADS] [CrossRef] [Google Scholar]
  32. Rosswog, S. 2015b, Liv. Rev. Comput. Astrophys., 1, 1 [CrossRef] [Google Scholar]
  33. Rosswog, S. 2020a, ApJ, 898, 60 [NASA ADS] [CrossRef] [Google Scholar]
  34. Rosswog, S. 2020b, MNRAS, 498, 4230 [NASA ADS] [CrossRef] [Google Scholar]
  35. Saitoh, T. R., & Makino, J. 2013, ApJ, 768, 44 [NASA ADS] [CrossRef] [Google Scholar]
  36. Springel, V. 2010, ARA&A, 48, 391 [NASA ADS] [CrossRef] [Google Scholar]
  37. Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649 [NASA ADS] [CrossRef] [Google Scholar]
  38. Tricco, T. S. 2019, MNRAS, 488, 5210 [NASA ADS] [CrossRef] [Google Scholar]
  39. Valdarnini, R. 2016, ApJ, 831, 103 [NASA ADS] [CrossRef] [Google Scholar]
  40. Wadsley, J. W., Keller, B. W., & Quinn, T. R. 2017, MNRAS, 471, 2357 [NASA ADS] [CrossRef] [Google Scholar]
  41. Wissing, R., & Shen, S. 2020, A&A, 638, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. 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:

(A.1)

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

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

(A.2)

The grad-h part of equation above can be estimated differentiating the constraint with respect ra which, after some algebra, gives:

(A.3)

when combined with expression (A.2) gives:

(A.4)

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

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

(A.5)

where Xb 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 Xb, as for example Xb = mb and Xb = 1 both reproducing the standard volume elements Vb = mb/ρb. Other choice is where P is the pressure and k ≤ 1. In this case, the estimator is strictly constant only in isobaric systems.

(A.6)

(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:

(A.8)

with

(A.9)

(b) Adaptive , being the standard SPH density. Because of the normalization condition, this choice leads to a constant kb = ∑cXcWbc(hb). Now we have

(A.10)

(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,

(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/∂ha). The result relies in the choice the estimator Xa used to compute the density.

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

    (A.13)

  • The choice requires a bit more algebra:

    (A.14)

    with:

    (A.15)

    and,

    (A.16)

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

    (A.17)

All Tables

Table 1.

Value of the different parameters in profiles A, B, and C that mimic different types of sharp density gradients.

Table 2.

L1 values for errors E1 and E2 at t = 1.5 in the hydrostatic square test.

Table 3.

Simulated models in the Kelvin–Helmholz test.

Table 4.

Simulated models in the wind-cloud test.

Table 5.

Summary of all calculated Sedov tests stating the values of σ, heat transport constant in the AV, VE choice, density peak, and L1 errors at t = 0.09.

Table 6.

Summary of the features of the approach presented in this work within the context of the works on which it is based.

All Figures

thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 4.

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

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 10.

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

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 14.

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

In the text
thumbnail Fig. 15.

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

In the text
thumbnail 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 (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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 18.

Sedov explosion. Evolution of the averaged L1 error, (Eq. (25)), in the shocked region corresponding to models S4,  (Xa = ma), and S8,  (), 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 (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.

In the text
thumbnail 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 .

In the text
thumbnail 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.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.