© The Authors 2026

1 Introduction

Smoothed particle hydrodynamics (SPH) is a Lagrangian particle method that has been applied to a wide range of topics within astrophysics. The traditional SPH method can be derived from the principle of least action and the Euler-Lagrange equations (Monaghan 1988; Price 2012), resulting in a numerical scheme for hydrodynamics that spatially conserves linear momentum, angular momentum, entropy, and energy exactly and that leaves the error in conservation mainly dependent on the time integration scheme¹. Many of the drawbacks of the traditional SPH method, such as handling shocks, shear flow, convergence issues, and density gradients, have been greatly improved over the years. In addition, there are still numerous aspects of the method that remain unexplored, with the potential for even more improvements in the future. In this article we focus on improving the convergence qualities of SPH.

Since SPH is a particle-based method, the fluid is discretized by mass elements rather than volume elements, as in grid-based codes. A smoothing kernel is used to interpolate the fluid quantities at any point in space. This includes the density, which determines the effective volume of the particles. While the sum of volume elements in a grid code remains constant, the sum of particle volumes in an SPH code is highly dependent on the accuracy of the interpolation. This accuracy is mainly determined by the smoothing kernel and how well the particles are distributed within the kernel. An optimally distributed set of particles satisfies, for each particle a, the following constraints: $$Q_{0,a}=\sum _b\frac{m_b}{{\rho }_b}{W}_{ab}=1$$(1) $$Q_{1,a}=\sum _b\frac{m_b}{{\rho }_b}({\mathbf{r}}_b-{\mathbf{r}}_a)W_{ab}=0,$$(2)

where the summation runs over neighboring particles, b, with mass m_b, density ρ_b, and kernel weighting function W_ab. Here Q₀ represent the partition of unity, and Q₁ indicates that there is no preferred direction or bias in the distribution of particles. The larger the deviation from these values, the worse the interpolation of fluid quantities. Related to this is the effective area between particles, which is analogous to the flux area between volume elements in grid-based codes. In a grid-based code, the areas fully enclose the volume of the resolution element (${\int }_s\hat{n}\mathrm{d}S=0$), where S is the bounding surface of the resolution element and $\hat{n}$ is the unit normal vector pointing outward from the surface. However, in SPH this is not necessary true and will again depend on the particle distribution such that there can be a non-canceling of particle areas (${\int }_s\hat{n}\mathrm{d}S\ne 0$). As momentum flux is made to be symmetric in SPH (such that momentum is conserved), non-canceling of particle areas leads to an automatic re-meshing mechanism in SPH. This causes particles to move to minimize this error and move toward a more optimal particle dis-tribution². This will occur even in the absence of hydrodynamic forces (constant pressure) and has thus come to be known as the zeroth order error of SPH. Particles that are optimally distributed fulfill the cancellation of areas and the following conditions: $$E_{0,a}^i=2\sum _b\frac{m_b}{{\rho }_b}\overline{\mathrm{\nabla }{W}_{ab}}=0$$(3) $$E_{1,a}^{ij}=2\sum _b\frac{m_b}{{\rho }_b}({\mathbf{r}}_b-{\mathbf{r}}_a{)}^i\overline{{\mathrm{\nabla }}_a^j{W}_{ab}}={\delta }^{ij}.$$(4)

Here E₀ is the zeroth order error³ and E₁ is the linear gradient error. The indices i, j ∈ {x, y, z} denote Cartesian vector components and δ^ij the Kronecker delta. Here we use the GDSPH gradient operator (Wadsley et al. 2017), which we continue to use for the rest of this paper. The momentum equation for SPH is derived from the Euler equations and given by $$\frac{\mathrm{d}v}{\mathrm{d}t}=\frac{-\mathrm{\nabla }P}{\rho }\to \frac{\mathrm{d}{v}_a}{\mathrm{d}t}=-\sum _b{m}_b\left(\frac{P_a+P_b}{{\rho }_a{\rho }_b}\right)\overline{\mathrm{\nabla }{W}_{ab}}.$$(5)

To get the error terms of this function, we can Taylor expand P_b around r_a : $$\begin{array}{rl}\frac{\mathrm{d}{v}_i}{\mathrm{d}t}& =-\frac{2P_a}{{\rho }_a}\sum _b\frac{m_b}{{\rho }_b}\overline{{\mathrm{\nabla }}_a^i{W}_{ab}}\\ & -\frac{2{\mathrm{\nabla }}_i{P}_a}{{\rho }_a}\sum _b\frac{m_b}{{\rho }_b}(r_b-r_a{)}^j\overline{{\mathrm{\nabla }}_a^j{W}_{ab}}+O(h^2).\end{array}$$(6)

It can be seen that our momentum equation is of the second order if the particle fulfills the E₀ and E₁ conditions. Most of the time this is not the case, but we can develop optimizations to try and minimize these kinds of errors. These optimization often involve different smoothing kernels, numbers of neighbors, corrections to the gradients, different gradient operators, or different choices of effective particle volume.

Linear-exact gradient corrections have become popular in recent times due to the increased accuracy they offer in subsonic flows and in shearing flows. Here, a matrix inversion of the E₁ condition is performed and then applied to the gradient of the kernel to ensure that the E₁ condition is fulfilled. There are a few variations to this correction. The most straightforward is to simply invert E₁ in Eq. (4) (Bonet 1999; Price & Monaghan 2004). In the CRKSPH formulation (Frontiere et al. 2017) and the SPH method presented by Rosswog (2025), the kernel and its derivatives are reproduced to enforce both the zeroth order and the linear order; however, when deriving the pairwise forces to uphold linear momentum conservation, the zeroth order correction is effectively lost, and only the linear correction remains. There might be benefits in these reproducing kernels over other linear correction methods, but that remains to be seen. Another way of constructing linear-exact gradients for SPH 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. Similarly to Bonet (1999); Price & Monaghan (2004), a matrix inversion is then performed to linearly correct for the gradients. This method is known as the integral smoothed particle hydrodynamics (ISPH) scheme (García-Senz et al. 2012). Integral-based estimates have been readily applied in the past for second derivatives, as they are much less noisy than analytical second derivatives of the smoothing kernel (Brookshaw 1985; Monaghan 2005). Moreover, ISPH have shown great advantages over regular SPH in subsonic flows (García-Senz et al. 2012; Rosswog 2015; Valdarnini 2016). The disadvantages of these linear corrections include the additional computational cost and the discrepancy in errors between the density estimate and velocity gradients, potentially leading to entropy errors. Global angular momentum is also no longer fully conserved in these formulations, as the pairwise force is not always radially aligned between particles. This becomes a resolution-dependent error, and it has been stated to remain relatively small (Frontiere et al. 2017); however, this error is dependent on both resolution and the given particle distribution. On the other hand, this is offset by more precise local angular momentum transport due to more accurate gradients.

The interpolation kernel is the foundation of the SPH method and is ever present in determining the accuracy of a simulation, even when applying the optimizations mentioned above. For a long time, B-splines were the most popular smoothing kernel to use due to their compact support, good interpolation, and flexibility in polynomial degree. However, an issue with these kernels has always been that they are susceptible to the pairinginstability. The pairing-instability occurs when the number of neighbors exceeds a certain critical value, causing particles to clump together and consequently reducing the resolution of the simulation (Thomas & Couchman 1992; Morris 1996; Børve et al. 2004; Price 2012; Dehnen & Aly 2012). This critical value will depend on the kernel, but for the cubic kernel, it lies in the region of N_smooth = 50. Increasing the number of neighbors leads to more accurate low-order errors (E₀, E ₁, Q₀, Q₁) and less particle noise. This is particularly desirable when modeling subsonic flows, where smooth velocity fields are essential⁴. This is why Wendland kernels have become very popular in recent times, as they are stable against the pairing instability at all neighbor numbers. It was shown in Dehnen & Aly (2012) that a positive Fourier transform is a necessary condition for stability against the pairing instability, though it is still unclear why negative values in the Fourier transform triggers the pairing instability. It has been hypothesized that this is a consequence of the particles trying to re-order themselves in order to minimize the total internal energy (U) given a fixed entropy. While the uniform particle distribution always represents a local minimum, if a paired particle distribution results in a lower internal energy, then the uniform particle distribution is only meta-stable. Cubic spline kernels are hypothesized to exhibit this meta-stable property, as it has an oscillating over-under estimation of density given a uniform particle distribution for varying neighbor numbers. The Wendland kernel and other kernels with positive definite Fourier transform seem to exhibit an overestimation of the density, which continuously decreases as number of neighbors increases. Dehnen & Aly (2012) hypothesized that pairing occurs when the following condition is fulfilled: ρ(N_smooth/f ) < p(N_smooth); for 1 < f ≤ 2. The overestimation of the density by the Wendland kernel can be especially high when using low neighbor numbers. Nevertheless, one can correct for the kernel bias by adjusting the self-interaction term (W(0, h)), thus effectively removing this bias from the density estimate while leaving the gradients unaffected. Though the Wendland kernels have been popularized, there are many other interpolation kernels that meet the criteria of having a positive Fourier transform, Gaussian shape, and compact support. These include the Wu family of kernels (Wu 1995), the missing Wendland kernels (Schaback 2011), and the Buhmann family of kernels (Buhmann 1998). Another family of kernels is the Sinc kernels, $W(q)={(\sin (\frac{\pi }{2}q)/(\frac{\pi }{2}q))}^n$ (Cabezón et al. 2008; García-Senz et al. 2014), which is directly linked to the Dirac-δ function. Sinc kernels, as with the B-spline family of kernels, can be subject to the pairing instability. However, the critical neighbor number for the instability can be pushed higher by increasing the sharpness of the kernel, which is controlled by the kernel exponent n. A consequence of increasing the sharpness is that it usually gives a worse interpolation at lower neighbor numbers. This is likely related to the larger difference in weight between the inner and outer regions produced by sharper kernels. Recently, a linear combination of two Sinc kernels was proposed by Cabezón & García-Senz (2024) to generate a kernel that is more resistant to the pairing instability while retaining good interpolation properties over a wide range of neighbors numbers. In this work, we aim to design smoothing kernels specifically optimized for SPH rather than using interpolation kernels developed for other purposes. Our approach involves optimizing a linear combination of kernels, similar to the method in Cabezón & García-Senz (2024), to produce the most accurate kernel for a given neighbor number. A big challenge in this work was determining the most natural initial conditions to minimize any kind of bias.

Lattice initial conditions are often used to provide an initially optimal distribution of particles. In this case, particles are put in a cubic or closed package lattice. At first this might seem to be a nice solution, as the E₀ and E₁ conditions are fulfilled. But there are many disadvantages to actually using a lattice for the initial condition. First, the lattice structure is easily broken by shocks and shear flows, which will generate aggressive noise as the particles move off the lattice structure, and this will cascade through the simulation volume. Second, particle lattice distributions can introduce directional bias and simulation artifacts. For instance, a shock wave propagating along a preferred direction of the lattice can disproportionately gather particles in that direction, leading to undesirable oscillations. Third, the two issues above are exacerbated when dealing with density gradients in the initial conditions. Lattice particle distributions can give both an artificial positive result for SPH, as it provides very high accuracy as long as the lattice holds up, but also an artificial negative result, as significant noise is generated when the lattice breaks down. A much more natural distribution for SPH is the glass distribution, which SPH will always strive to re-mesh toward as particles are disturbed by such forces as shocks and shear flows. Glass distributions are often obtained by letting the particles relax, beginning with a lattice or random distribution, adding some random velocities, and using a velocity damping term to eventually approach a steady state. This can be done with the influence of gravity to generate collapsed density structures or without it to generate uniform distributions. However, this method has limitations and downsides. First, one cannot generate density contrasts for setups that are purely hydrodynamic⁵ (the Kelvin-Helmholtz setup for example). Second, glass generation performed in this way can be costly, as oscillations can take a long time to damp. Third, damping the velocities too quickly can lead to inaccurate distributions. In this paper, we have instead developed a glass generator that relies on the SPH method to find a relaxed glass distribution for arbitrary density gradients.

The use of lattice initial conditions is especially problematic when it comes to assessing the quality of a smoothing kernel, as the results do not reflect the natural distribution of particles and give a positive bias toward smoother kernels, which usually have larger E₀ errors in glass distributions (given the same neighbor count). In Dehnen & Aly (2012), the correction for the kernel bias was made based on the bias when interpolating a closed-packed lattice distribution. The bias was also given by a simple power law. In this paper we calculate a new correction to the kernel bias for all our kernels that is based on the generated glass that each kernel itself produces. We also give a piecewise function to capture the kernel bias more accurately over a wider-range of neighbor numbers.

In Sect. 2 we present all the smoothing kernels used in this paper. In Sect. 3 we go through the methods used in this paper, including the SPH methods, the IC generator, adjustment of density bias, and the optimization strategy. In Sect. 4 we present the results from the tests. And finally in Sect. 5 we discuss our results.

2 Kernel properties

All the kernels we go through in this section go to zero at $$q=(r/2h)>1,$$(7)

even if not stated explicitly in the equations. All kernels are normalized with $$\sigma =\frac{1}{{\int }_0^1{\int }_0^{\pi }{\int }_0^{2\pi }W(r)r^2\sin (\theta )\mathrm{d}r\mathrm{d}\theta \mathrm{d}\psi }.$$(8)

A useful tool for analyzing kernels is their Fourier transform. Despite having similar shapes, kernels can exhibit quite different Fourier transforms. As previously mentioned, negative values in the Fourier transform are responsible for the pairing instability. The lower the wave numbers at which these negative values appear, the lower the critical number of neighbors required to trigger this instability. The Fourier transform is defined as $$F_3[W(r)](k)=4\pi k^{-1}{\int }_0^{\mathrm{\infty }}\sin (kr)W(r)r\mathrm{d}r,$$(9)

where W(r) is the smoothing kernel and k is the wave number. All the kernels that we mention in this section are positivedefinite, which means that they have as strictly positive Fourier transform. Positive definite kernels all produce a general overbias in their density estimation, which become more prominent at low neighbor numbers, as the bias is dominated by the self contribution. The neighbor number at which this bias becomes negligible depends on the smoothing kernel (higher order/more peaked kernels give higher bias). This self contribution can be corrected for as described in Sect. 3.3.

2.1 Generalized Wendland kernel

The Wendland kernels is a family of compactly supported radial functions generated by a dimension walk (integration) of a truncated power function (Wendland 1995). $$\begin{array}{r}{\psi }_l(r)=(1-r{)}_{+}^l,\\ l(k)=\lfloor 1.5+k\rfloor +1,\\ {\varphi }_k(r)={\mathcal{I}}^k{\psi }_l.\end{array}$$(10)

Here k is the smoothness parameter, l(k) is the order of the kernel, and $\lfloor x\rfloor$ is the floor operator. We set l(k) so that the generated function is strictly positive definite and radial on R³ (following D ≤ 2l - 1 for D dimensions). The Wendland kernels produce the minimal polynomial degree for a given space dimension and smoothness. While initially only determined for positive integer smoothness parameters, this property has been generalized for non-integer smoothness parameters in the generalized Wendland functions (Chernih & Hubbert 2014): $$\begin{array}{r}{\varphi }_k(r)=\frac{\mathrm{\Gamma }[l(k)+1]}{2^{l(k)+k}\mathrm{\Gamma }[l(k)+k+1]}{(1-r^2)}^{l(k)+k}{r}^{-l(k)}\\ \times {}_2{F}_1(\frac{l(k)}{2},k+\frac{l(k)+1}{2},l(k)+k+1,1-\frac{1}{r^2}).\end{array}$$(11)

Here Γ is the gamma function and ₂F₁ is the hyper-geometric function. The family of Wendland kernel functions for 3D is then defined as $$W(r)=\sigma {\varphi }_k(r).$$(12)

Positive integers of k represent the classic Wendland kernels, where the resulting kernels have polynomial structure. The Wendland kernels with positive integers used in this paper (k = 0,1,2,3) are $$W_{C0}(q)=\sigma (1-q{)}^2,$$(13) $$W_{C2}(q)=\sigma ((1-q{)}^4(1+4(1-q))),$$(14) $$W_{C4}(q)=\sigma ((1-q{)}^6(1+6(1-q)+\frac{35}{3}(1-q{)}^2)),$$(15) $$W_{C6}(q)=\sigma ((1-q{)}^8(1+8(1-q)+25(1-q{)}^2+32(1-q{)}^3)).$$(16)

Positive half-integers are known as the missing Wendland kernels (Schaback 2011). They are compactly supported and polynomial but with additional logarithmic and square-root terms ( k = 0.5): $$\begin{array}{r}{W}_{C{M}_{05}}=\sigma (\sqrt{1-q^2}+6.5q^2\sqrt{1-q^2}\\ +1.5q^2(4+q^2)(\log q-\log (1+\sqrt{1-q^2}))).\end{array}$$(17)

Positive non-integer kernels can also be generated; these involve a more complicated form (k = 0.4): $$W_{C{M}_{04}}=\sigma \frac{0.375(1-q^2{)}^{2.4}{}_2{F}_1(1,1.9,3.4,1-\frac{1}{q^2})}{q^2}.$$(18)

To avoid handling hyperbolic functions within our code we approximate this function with a polynomial fit. Higher-order kernels have smoother derivatives, which in combination with higher neighbor numbers act to decrease the sensitivity to particle disorder. The trade-off of using a higher-order kernel is that it becomes less accurate than its lower-order counterpart at low neighbor numbers (Monaghan 1992; Rosswog 2015). We refer to the Wendland kernels as C2, C4, C6, CM05, and CM04 throughout this paper.

2.2 Wu kernel

The Wu approach starts with the function (Wu 1995) $$\varphi =(1-r^2{)}_{+}^l.$$(19)

Another function is then generated by convolution: $${\varphi }_{\ell }(r)={\int }_{-1}^1(1-t^2{)}^{\ell }(1-(2r-t{)}^2{)}^{\ell }\mathrm{d}t.$$(20)

Using this function, the Wu family of kernels is generated by dimension walk (derivative): $${\varphi }_{l,k}={\mathcal{D}}^k{\varphi }_l.$$(21)

The Wu kernel (φ_l,k) is strictly positive definite and radial in R^d for D ≤ 2k + 1, with polynomial of degree 41 - 2k + 1 and smoothness of C^2(l-k), which indicates how many times the kernel is continuously differentiable (2(l - k) times). The polynomial degree is higher for the Wu kernel than the Wendland kernel at a certain smoothness. In this paper we have chosen to look at the φ_1,2 and φ_1,3: $$W_{WU2}=(1-q{)}_{+}^4(4+16q+12q^2+3q^3),$$(22) $$W_{WU4}=(1-q{)}_{+}^6(6+36q+82q^2+72q^3+30q^4+5q^5).$$(23)

For a prescribed smoothness the polynomial degree of Wendland kernels are lower than that of Wu kernels. An interesting thing about Wu kernels is that their Fourier transform is positive, decays quickly, and goes to zero in periodic intervals. We refer to the Wu kernels as WU2 and WU4 throughout this paper.

2.3 Buhmann kernel

Buhmann’s approach (Buhmann 1998) involves integrating a positive function f(t) = t^α(1 - t^δ)^ρ₊ against a strictly positive definite kernel $K(t,r)={(1-\frac{r^2}{t})}_{+}^{\lambda }$. Buhmann’s family of kernels thus take the general form $$\varphi (r)={\int }_0^{\mathrm{\infty }}(1-r^2/t{)}_{+}^{\lambda }{t}^{\alpha }(1-t^{\delta }{)}_{+}^{\rho }\mathrm{d}t.$$(24)

Here 0 < δ ≤ 0.5, ρ > 0.5 and the kernel is strictly positive and radial in R^D ≤ 3 for λ ≥ 0 and $-1<\alpha <\frac{\lambda -1}{2}$. Buhmann’s family of kernels consists of a polynomial and an additional logarithmic term. The Buhmann kernel that we include in this paper is (α = $\delta =\frac{1}{2}$ ρ = 1 λ = 2 from Buhmann 2003): $$W_{Buh}(r)=\sigma (12q^4\log q-21q^4+32q^3-12q^2+1).$$(25)

We refer to this kernel as BUH throughout this paper.

3 Method

3.1 SPH methods

We performed the optimization using both the GDSPH and ISPH methods. We did this because the “best” kernel may differ between the two methods due to the removal of the linear error term. The density of both GDSPH and ISPH was calculated using the standard density estimate $${\rho }_a=\sum _b{m}_b{W}_{ab}(h_a).$$(26)

We used the same artificial viscosity scheme for all the simulations. To keep things as simple as possible, we implemented a classic SPH artificial viscosity (Monaghan 1992) with constant α and β. We used α = 1 and β = 2 for the shock cases (Sod-shocktube and Evrard collapse), and for the subsonic cases (Gresho-Chan, Kelvin-Helmholtz, and the subsonic blob test) we used α = 0.05 and β = 2.

3.1.1 GDSPH gradient

We implemented the GDSPH gradient operators of Gasoline2 (Wadsley et al. 2017) in this paper. This variant of SPH is identical to regular SPH when the density is uniform, but it improves gradient accuracy in cases involving sharp density gradients. The momentum equation and internal energy equation in GDSPH is written as: $$\frac{\mathrm{d}{\mathbf{v}}_a}{\mathrm{d}t}=-\sum _b{m}_b\left[\frac{f_a{P}_a+f_b{P}_b}{{\rho }_a{\rho }_b}\overline{{\mathrm{\nabla }}_a{W}_{ab}}\right],$$(27) $$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\frac{f_a{P}_a}{{\rho }_a}\sum _b\frac{m_b}{{\rho }_b}{\mathbf{v}}_{ab}\cdot \overline{{\mathrm{\nabla }}_a{W}_{ab}}.$$(28)

Here, f_a are entropy correction terms to remove the linear entropy error caused by a mismatch between GDSPH operators and the density estimate (see Wadsley et al. 2017).

3.1.2 Integral based gradient

In ISPH (García-Senz et al. 2012), the gradients are calculated from an integral expression so that there is no need to explicitly calculate the analytic derivative of the kernel. The kernel gradients are also linear-corrected by the use of a matrix inversion of the E₁ error. Effectively, we replace the kernel gradients of Eqs. (27) and (28) with ($\overline{{\mathrm{\nabla }}_a{W}_{ab}}\to \overline{G_{ab}^k}$): $$(G_{a,j}{)}^k=\sum _{d=1}^3{C}_a^{kd}{r}_{ba}^d{W}_{ab}(h_j),$$(29) $$(G_{b,j}{)}^k=\sum _{d=1}^3{C}_b^{kd}{r}_{ba}^d{W}_{ab}(h_j),$$(30)

where j is either a or b. $$\overline{G_a^k}=0.5(G_{a,a}^k+G_{a,b}^k),$$(31) $$\overline{G_b^k}=0.5(G_{b,a}^k+G_{b,b}^k).$$(32)

The linear correction matrix, C, is calculated with $$C_a^{ki}={\left(\sum _b\frac{m_b}{{\rho }_b}{r}_{ba}^k{r}_{ba}^i\overline{{\mathrm{\nabla }}_a{W}_{ab}}\right)}^{-1}.$$(33)

This form is slightly different than the one proposed by previous authors(García-Senz et al. 2012; Rosswog 2020). This is because it involves the symmetric kernel $\overline{{\mathrm{\nabla }}_a{W}_{ab}}$ within the linear correction matrix, instead of the one-sided kernel gradient ${\mathrm{\nabla }}_a{W}_{ab}(h_a)$. The momentum and internal energy equations become $$\frac{\mathrm{d}{\mathbf{v}}_a}{\mathrm{d}t}=-\sum _b{m}_b\left[\frac{P_a\overline{G_a^k}+P_b\overline{G_b^k}}{{\rho }_a{\rho }_b}\right],$$(34) $$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\frac{P_a}{{\rho }_a}\sum _b\frac{m_b}{{\rho }_b}{\mathbf{v}}_{ab}\cdot \overline{G_a^k}.$$(35)

This matrix inversion requires us to do an extra neighbor loop and requires us to save six additional quantities per particle (symmetric matrix). This linear correction can introduce a global angular momentum error as pairwise torques can be generated between particles. This is because the forces between particles are no longer ensured to be aligned along the line connecting particles. This is in general true for any kernel correction of order greater than zero. This effectively results in a resolution/kernel dependent error in global angular momentum conservation.

3.2 Initial condition: Glass generator

As discussed in the introduction, producing natural and nonbiased initial conditions are important for SPH. Of particular difficulty is the generation of non-uniform density distributions. Lattice initial conditions can be used in non-uniform density distributions, by either varying particle mass (Lombardi et al. 1999) or stretching a uniform lattice (Rosswog 2009; Price et al. 2018). The main issue with these methods are the use of the lattice initial condition, which as discussed are artificial, biased and experience excessive noise when the lattice breaks. Varying particle masses are undesirable as well, due to particles having more difficulty in balancing the particle distribution and generally leading to larger gradient errors and noise. In Diehl et al. (2015) a more general initial condition generator is proposed, using the concept of weighted Voronoi tessellation. This works well in generating non-uniform density distributions. The main issue with this approach is that the generated distribution does not necessarily correspond to the same natural particle distribution that a given SPH method would produce, as it does not take into account the smoothing kernel or the gradient operators. This means that it will also represent a slightly biased initial condition. In Rosswog (2020), a more SPH-like version of this method is applied, where the SPH momentum equation is used, together with an artificial pressure based on the current density error. This is an excellent approach for generating non-biased initial conditions for SPH, as it ensures a near noise free distribution at the start of your simulation (as one is using the same momentum equation to evolve the system). We propose a similar method in this paper, but add some important alterations to ensure that both local and global errors are minimized. In addition, this method will relax to a given density profile, no matter what your starting initial condition is (given that the mass of the system is correct). Similar to Rosswog (2020), we used an iterative method where during each iteration, i, we updated the position directly for each particle: $${\mathbf{r}}_{a,i+1}={\mathbf{r}}_{a,i}+\mathrm{\Delta }{\mathbf{r}}_{a,i}$$(36)

The displacement of each particle is determined by the following equations: $$\mathrm{\Delta }{r}_{a,i}=K_{global,i}\frac{\mathrm{d}{v}_a/\mathrm{d}t}{f_{max,i}}{h}_a$$(37) $$f_{max,i+1}=max\left(\left|\frac{\mathrm{d}v}{\mathrm{d}}\right|\right).$$(38)

The term K_global,i will determine the largest multiple of the smoothing length that a particle can be displaced. The variables $\frac{\mathrm{d}v}{\mathrm{d}t}$ and ρ_i are calculated with the chosen density estimate and gradient operator (in this paper GDSPH (Eq. (27)) and ISPH (Eq. (34))). Hence the produced glass will represent the relaxed distribution of the chosen gradient operator, smoothing kernel and neighbor number. The pressure in the momentum equation is replaced by an artificial pressure given by $$P_a={\left(\frac{{\rho }_{tar}}{{\rho }_a}\right)}^n.$$(39)

Here ρ_tar is the target density at the current position and ρ_a is the estimated density for the particle (Eq. (26)). Higher n means that fitting to the density is more important than balancing the particle distribution, that is, reducing E₀. This is because a minimization in density error does not necessary mean a minimization in E₀. As such, density errors and particle distribution errors will drive the particles into a position that minimizes these errors, just as the case for SPH naturally⁶. The term K_global,iis adjusted after each step following the instructions outlined below:

• If the step leads to an overall worse mean $\sqrt{\sum {\left(\frac{\mathrm{d}{v}_i}{\mathrm{d}t}\right)}^2}$, then decrease $K_{global,i+1}=\frac{1}{d_1}{K}_{global,i}$.

• If the step leads to an overall better mean $\sqrt{\sum {\left(\frac{\mathrm{d}{v}_i}{\mathrm{d}t}\right)}^2}$, then increase $K_{global,i+1}=d_2{K}_{global,i}$.

• Relaxation requires that d₁ > d₂ > 1.0. We set d₁ = 1.8 and d₂ = 1.1 during the global relaxation and d₁ = 1.1 and d₂ = 1.002 during local relaxation.

• The global relaxation stage quickly reduces the step size to interparticle size to get a good gradient estimate of the current global density profile. The step size is then increased again to a factor of the initial (C_loop,nK_global,0). Here we set C_loop,0 to 0.5 and reduce it by half each time C_loop,n = 0.5C_loop,0 if there is no improvement to the global density error. If one starts from an IC that captures the global IC well, then this step is not necessary, and only local relaxation is needed.

• In the local relaxation stage, which occurs when global relaxation stage is finished, the K_global,i in Eq. (37) is near the average interparticle scale. From here on we slow down the rate of relaxation to ensure a more accurate local particle distribution. We stop the IC generator when max(∆r_a/h_a) < 2 ∙ 10⁻³, at which point particle motions are negligible, and we consider the distribution relaxed.

The initial value of K_global,i is a parameter that can be chosen, but by default it is set to half the system size, as this ensures that the desired density profile is obtained regardless of the input IC. A much lower value can be used if one has an IC that is close to the desired density profile.

Results from the IC generator in the case of the ∆ρ = 8, n_x = 64, N_smooth = 64 Kelvin-Helmholtz IC of Sect. 4.3 can be seen in Fig. 1. From this we can see that we have a proper glass structure in the particle distribution that follows the analytical density profile very closely. As different relaxed glasses will be produced by different smoothing kernels and neighbor numbers, we run the relaxation process for each new kernel and neighbor number when generating ICs in this paper.

3.3 Adjusting density bias

Before we optimized our linear combined kernel, we adjusted the self-bias of all our kernels. This mainly effects the density, giving us a corrected density estimate: $${\rho }_a={\rho }_{a,estimate}-ϵm_{a}W(0,h_a).$$(40)

This density correction is a constant and is applied uniformly to all particles. We adjust e so that the density of a relaxed hydrostatic glass is as close as possible to ρ = 1.0 (given that M_tot/V_tot = ρ = 1.0) for a wide range of neighbor numbers. A uniform glass produced by the IC generator will locally be independent of the resolution and only be a function of the number of neighbors and the given kernel. As such we can run low resolution cases with our IC generator and adjust the kernel bias until we produce a glass with a mean density error of around ρ_err ≈ 10⁻⁵, which is much less than the general SPH density noise. We run this for a wide range of neighbor numbers (N_smooth ∈ {16,24,32,48,64,96,128,256,512}). After 512 neighbors, the kernel bias is set to 1. We could then interpolate between ranges using a linear piecewise function to correct the kernel: $$ϵ(N_{smooth})=\frac{c_{i+1}-c_i}{N_{min,i+1}-N_{min,i}}(N_{smooth}-N_{min,i})+c_i.$$(41)

Here, c_i is the e value at N_smooth = N_min,i. This was performed for both GDSPH and ISPH. The coefficients and neighbor number can be found in Table B.1 for GDSPH and Table B.2 for ISPH.

Fig. 1 Relaxed glass IC generated by the IC generator outlined in Sect. 3.2. Here we model the IC of the Kelvin Helmholtz instability (Sect. 4.3), with a density contrast of ∆ρ = 8. Top : particle distribution within a thin slice of thickness ∆ z = 0.015. Bottom : density profile in the y direction. The red line depicts the analytical density profile and the black dots the density of the particles. The glass generated density profile is in excellent agreement with the analytical solution.

3.4 Optimization strategy

In this method, we generate a new kernel by linear combination: $$\begin{array}{rl}{W}_{new,i}({\alpha }_i,{\alpha }_{i+1})& =({\alpha }_i{W}_{new,i-1}+(1-{\alpha }_i)W_i){\alpha }_{i+1}\\ & +(1-{\alpha }_{i+1})W_{i+1}.\end{array}$$(42)

Here, α_i and α_{i + 1} are the coefficients that we optimized for. The process starts with i = 1 and an initial kernel W_new,0. After determining α_i, α_{i + 1} by minimizing some cost function, we generated a new kernel (W_new,i). We then repeated this process with the newly generated kernel and two other kernels. This iterative process continued until we used all the kernels that we wished to include in our linear combination. We used the differential evolution method to find the global minimum of the two parameters for each iteration. We set the initial parameter value to (α_i,α_i+1) = [0.0,0.0], and we set the minimum and maximum parameter boundaries to be bnds = [(−3.0, 3.0), (−3.0,3.0)].

After the optimization was completed, we fit the linearcombination of kernels to a polynomial in the form $$P(x)=\sum _{i=0}^8{c}_i{x}^i$$(43) $$W(q)=\sigma (1+q^{2}P(q)).$$(44)

During the work of this paper we used two different cost functions and two different cases to optimize the kernel for. These are described in the sections below.

3.4.1 Optimized for hydrostatic glass

Perhaps, one of the simplest cases to determine the quality of a smoothing kernel comes from it’s ability to generate a hydrostatic glass that minimizes the leading errors of SPH. The hydrostatic glass models a homogeneous system in pressure equilibrium. A cubic box with length L = 1 and total mass M_box = 1 is setup. Particles are randomly distributed throughout the box. Due to the inherent re-meshing property of SPH, the particles will be reorganized to a glass-like structure by the low order errors (Eqs. (3) and (4)). At some point the distribution will be organized into a relaxed glass, where the errors have been minimized (there is minimal/oscillatory change in the mean errors from this point onward). The distribution and the magnitude of the errors in this final relaxed state will depend on both the smoothing kernel, number of neighbors and the SPH method. The leading order error (E₀) will locally be independent of the resolution (for uniform density), so we use a relatively low resolution of N = 32³ particles. We hypothesized that the relaxed glass with the smallest errors/noise would provide the best result in dynamic cases as well. We defined our cost function for this case as $$F=\delta {\rho }_{err}+E_0+0.25E_1.$$(45)

Here E₀ and E₁ are given by the norm of the error vectors defined in (Eqs. (3) and (4)) respectively. The term δρ_err is given by $$\delta {\rho }_{err}=\frac{1}{N}\sum _i|{\rho }_i-{\rho }_{mean}|$$(46)

Each kernel has had its self-bias adjusted, so that ρ_mean will be within 10⁻⁵ of ρ = 1.0, for their respective glass and kernel (see Sect. 3.3), such that we accurately measure the density noise. We only optimized this case with the GDSPH formula, and we refer to these kernels as O_{glass,Nsmooth} (N_smooth = 32, 64, 128, 256) for the rest of the paper. The resulting polynomial coefficients for the optimization are given in Table C.1.

The kernel function, its derivative and its Fourier transform can be seen in Fig. 2 for the optimized hydrostatic glass kernels together with WU₂, BUH, CM₀₄, CM₀₅ kernels. All the optimized kernels for different N_smooth look very similar to each other, and they are a bit steeper than the WU₂ kernels. However, some differences can be seen in the Fourier transform, where O_glass,32 has quite a flat trend (similar to BUH kernel), the higher N_smooth kernels experience stronger oscillations (behavior somewhere in-between the WU₂ kernel and the BUH kernel). Even though we allow for the linear combination to result in negative Fourier transforms, we only end up with one kernel that has negative Fourier transform in the high wave number regime (O_glass,64). The performance of these kernels are discussed more in the next section.

Fig. 2 Top : optimized kernel functions for a glass distribution (Sect. 3.4.1) together with the missing Wendland kernels (CM04, CM05), the Buhmann kernel (BUH), and the Wu kernel (WU2). Middle : corresponding derivative of these kernel function. Bottom: fourier transform of these kernels, scaled by k/10. The optimized kernels are quite similar, with some small variations being seen in the Fourier transform. The Oglass,32 kernel deviates the most from the other kernel, with less deep troughs in the Fourier transform. These kernel’s generate excellent glass distributions but degrade significantly in dynamical simulations.

3.4.2 Optimized for Gresho-Chan vortex

Initially, we assumed that a smoothing kernel with minimal errors in its final relaxed glass state would also yield the best results in dynamic simulations, but this turned out to not be the case. As different kernels have different sensitivity when it’s glass gets disturbed. This was clearly seen when testing out the Gresho-Chan vortex test case. Here, even though the density is uniform, the particle distribution is continuously disturbed by a shear flow.

The Gresho-Chan vortex simulates an inviscid fluid vortex in force equilibrium (Gresho & Chan 1990). The setup is as follows: We setup a 3D thin periodic box with length $L=(1,1,2\frac{\sqrt{6}}{n_x})$, here n_x is the average number of particles in the x direction, and the z boundary is set to be roughly 24 particle spacings. The pressure profile is setup following $$P(r)=\{\begin{array}{ll}5+\frac{25}{2}{r}^2& \text{for\hspace{0.17em}}0\le r<0.2\\ 9+\frac{25}{2}{r}^2-20r+4\ln 5r& \text{for\hspace{0.17em}}0.2\le r<0.4\\ 3+4\ln r& \text{for\hspace{0.17em}}r\ge 0.4\end{array},$$(47)

and the velocity profile follows $$v_{\theta }(r)=\{\begin{array}{ll}5r& \text{for\hspace{0.17em}}0\le r<0.2\\ 2-5r& \text{for\hspace{0.17em}}0.2\le r<0.4\\ 0& \text{for\hspace{0.17em}}r\ge 0.4\end{array}.$$(48)

Even though conservation of angular momentum within SPH is near exact, the test is challenging for SPH due to the artificial viscosity induced by the particle noise and the errors in the linear gradient (Eq. (4)). The artificial viscosity will cause angular momentum transport, which will lead to mass accretion toward the center and a loss in kinetic energy. Increasing the number of neighbors (Read & Hayfield 2012; Rosswog 2015), using artificial viscosity switches (Cullen & Dehnen 2010), and doing linear gradient corrections (García-Senz et al. 2012; Valdarnini 2016; Cabezón et al. 2017) all improve the convergence of this case.

While testing the Gresho-Chan vortex we found that only O_glass,32 outperformed the other kernels for the N_smooth it had been optimized for, while all the other O_glass kernels performed rather badly (can be seen in Fig. 5). Similarly, we saw bad result for all kernels that exhibit strong oscillatory behavior in their Fourier transform, such as the WU₂ kernel (Fig. 2). Therefore we decided to additionally optimize kernels based on results from the Gresho-Chan vortex following the procedure laid out in the previous sections. The main difference is the cost function, which we set to be the L_1,error at t = 1: $$F=L_{1,error}(t=1),$$(49) $$L_{1,error}=\frac{1}{N}\sum _{i=1}^N|v_{\varphi ,i}^a(r)-v_{\varphi ,i}^d(r)|$$(50)

Here $v_{\varphi ,i}^a$ stands for the analytical, and $v_{\varphi ,i}^d$ is the result from the simulation. Similarly to all other setups, we used the IC generator (Sect. 3.2) to generate an initial glass for each kernel and neighbor number before each run. We optimize both GDSPH and ISPH for this case; here the kernels optimized with GDSPH are given by O_Nsmooth and with ISPH O_I,Nsmooth, which was done for N_smooth = 32, 64, 128, 256. The kernels and the resulting coefficients that were used for the optimization are given in Tables C.2 and C.3.

The kernel function, its derivative and its Fourier transform can be seen in Fig. 3 for the optimized kernels together with C₂ and C₄ kernels. All the resulting kernels have similar amplitude and form as the Wendland C2 kernel (except for O_I,32), but with a minima in its derivative being a bit offset from the C2 kernel. Looking at the Fourier transform we can see quite a clear trend with both the O_Nsmooth and O_I,Nsmooth as we increase N_smooth ∙ Both O₃₂ and O_I,32 drop the fastest within k ∝ 0 - 3π but remain higher than the other kernels as we go to larger wave numbers. As we increase N_smooth we have a higher power between k ∝ 0 - 3π but lower for higher wave numbers. We have two resulting kernels with regions of negative Fourier transform. The O_I,32 kernel, which goes negative in the region k ∝ 4 - 5π, and becomes susceptible to the pairing instability as we go to higher N_smooth (we see instability at N_smooth = 128 and N_smooth = 256). The O₂₅₆ also has regions of negative Fourier transform at much higher wave numbers (k ∝ 10 - 12π). However, we did not see any sort of instability occurring in this kernel, at least up to N_smooth = 512. The performances of these kernels are discussed more in Sect. 4.

4 Results

In the following section we discuss the performance of these optimized smoothing kernels on a select number of test cases. First we go through the hydrostatic glass and the Gresho-Chan vortex, both of which we outlined and used for optimization of the kernels in Sects. 3.4.1 and 3.4.2. We then go through four additional tests to see how well the optimized kernels perform in general by testing the Kelvin-Helmholtz instability for two different density contrasts (∆ρ = 2, ∆ρ = 8), the classic Sod shocktube test, Evrard collapse and the subsonic blob test. The hydrostatic glass and Gresho-Chan vortex are run for all the optimized kernels. The ones optimized to hydrostatic glass with GDSPH: O_glass,32, O_glass,64, O_glass,128, O_glass,256. The ones optimized for Gresho-Chan vortex with GDSPH O₃₂, O₆₄, O₁₂₈, O₂₅₆ and with ISPH O_I,32, O_I,64, O_I,128, O_I,256. We also compare with a select number of kernel’s mentioned in Sect. 2 (CM₀₄, WU₂, BUH, CM₀₅, C2, C4, C6). This is done for N_smooth = 32,64,128,256. For the Kelvin-Helmholtz instability, Sod shocktube, Evrard collapse and subsonic blob test we compare the best-performing kernel at a given N_smooth with the popular Wendland C2 and C4 kernel, for both GDSPH and ISPH. The optimized kernels was implemented in both the Gasoline2 code (Wadsley et al. 2017) and the ChaNGa code (Menon et al. 2015). Results are from simulations performed with Gasoline2, but we have also confirmed that ChaNGa gives the same results in many of these cases.

4.1 Hydrostatic glass

The setup of this test case is outlined in Sect. 3.4.1. The system is evolved for a long enough time to reach a relaxed glass state. The quality of the resulting glass is quantified with the SPH errors (E₀ Eq. (3) E₁ Eq. (4)) and the density error (Eq. (46)). These are given for all kernels and chosen neighbor numbers in Fig. 4. While the kernel’s optimized for the hydrostatic glass perform in average the best considering the cost function (Eq. (45)), we can see that it loses out to some other kernels for specific errors. All the GDSPH optimized kernels for the Gresho-Chan vortex perform fairly well, similar to the BUH, CM₀₅ and C₂ kernel for the δρ_err and E₀ error. The O kernels does in general have smaller E₁ errors than the standard kernels. The only exception is the O₃₂ kernel that has higher E₁ error but has the smallest density and E₀ error out of all the kernels for N₃₂. The WU₂ kernel performs very well for higher neighbor numbers and is close to the results of the optimized kernel O_glass,128 and O_glass,256. The biggest improvement that can be seen when comparing the O_glass and O kernels with the Wendland C2 kernels, is in the E₁ errors that are significantly reduced. We can see that the O_I,₃₂ kernel performs very badly when relaxing to a hydrostatic glass with the GDSPH gradient operator, even at N_smooth = 32. This is, however, not the case when relaxing with ISPH at N_smooth = 32, indicating that glass relaxation is different when using ISPH.

Fig. 3 Top : optimized kernel functions for the Gresho-Chan vortex (Sect. 3.4.2), for both GDSPH and ISPH together with the popular Wendland kernels (C2, C4). Middle : corresponding derivative of these kernel function. Bottom : fourier transform of these kernels, scaled by k/10. The optimized kernel exhibit a clear trend in it’s Fourier transform, where the kernels optimized for higher Nsmooth emphasizes lower power beyond k > 3π. The general form of the optimized kernel is quite similar to the Wendland C2 kernel. In general GDSPH produces lower power than ISPH in the high wave number regime for a given Nsmooth. The OI,32 kernel is the only optimized kernel that exhibits a negative Fourier transform at low wavenumbers. Among the other optimized kernels, only O256 has a negative Fourier transform, but this occurs at higher wavenumbers and remains stable to pairing instability up to at least Nsmooth < 500.

Fig. 4 Relative increase or decrease of the δρerr (top), E0 (middle), and E1 (bottom) of the hydrostatic glass when relaxed for different kernels in respect to the corresponding error of the C2 kernel. The simulations here were performed with only the GDSPH method. The terms N32, N64, N128, and N256 refer to the number of neighbors: Nsmooth = 32,64,128,256, respectively. The δρerr, E0 and E1 of the C2 kernel for each Nsmooth is δρerr = 5.539 × 10−3,3.747 × 10−3,1.243 × 10−3,5.860 × 10−4, E0 = 2.763 × 10−3,2.024 × 10−3, 6.797 × 10−4,4.253 × 10−4, and E1 = 2.210 × 10−2,1.484 × 10−2,2.557 × 10−3,8.305 × 10−4. The optimized kernels can in general be seen to drastically reduce the E1 error relative to the C2 kernel. For the glass optimized kernel’s we can also see a general decrease in both δρ and E0 for their given Nsmooth.

Fig. 5 Relative increase or decrease of the L1,error of the Gresho-chan vortex at t = 1 and nx = 64 for different kernels in respect to the L1,error of the C2 kernel. The terms N32, N64, N128, and N256 refer to the number of neighbors: Nsmooth = 32,64,128,256. Left : simulations done with the GDSPH method. Right : simulations done with the ISPH method. The L1,error of the C2 kernel for each Nsmooth is L1,error = 1.49375 × 10−1, 9.57342 × 10−2, 4.70888 × 10−2, 2.99816 × 10−2 for GDSPH and L1,error = 9.34635 × 10−2, 4.66192 × 10−2, 2.64744 × 10−2, 2.22268 × 10−2 for ISPH. An improvement for the optimized kernel at a given Nsmooth across all Nsmooth can be seen for both methods, where the biggest improvements are for the N64 and N128 for GDSPH and N32 and N64 for ISPH. Of the non-optimized kernels, we observed that the C2 kernel performs best, except for Nsmooth = 32 (GDSPH), Nsmooth = 32,64 (ISPH), where CM05 performs better and C4 performs slightly better at Nsmooth = 256 (ISPH).

4.2 Gresho-Chan vortex

The setup of this test case is outlined in Sect. 3.4.1. We run the Gresho-Chan vortex up to t = 3. We also perform a convergence study, to see if the optimized kernels still perform well as we increase/decrease the resolution, as the optimization was performed at n_x = 64 resolution. For this convergence study we only compare with the Wendland C₂ and C₄ kernels. The optimization of the ON smooth kernels was done at t = 1, we can see the result of the L_1,error (see Eq. (50)) at t = 1 in Fig. 5 for both GDSPH and ISPH. From this table we can see that the optimized kernels for the Gresho-Chan vortex perform the best for their given N_smooth. The difference in L_1,error is more pronounced for lower N_smooth and for the GDSPH method compared to the ISPH method. A comparison of the azimuthal velocity structure (n_x = 64, t = 3) between the optimized kernel and the Wendland C₂ and C₄ kernel can be be seen in Fig. 6 for GDSPH and ISPH. For GDSPH we can see that the optimized kernel mitigates the decay of the vortex structure (in particular for N_smooth = 64). For ISPH the effect of the optimized kernels can mainly be seen as less noise in the vortex structure, while the general decay is similar across all the kernels (as this is highly dependent on the linear gradient estimate). A convergence study at t = 3 of the L_1,error can be seen in Fig. 7 for GDSPH and Fig. 8 for ISPH. Looking at the convergence study, we can see that the improvement of the optimized kernels in general stays true for higher and lower resolution.

Looking at the GDSPH results, we can see that while the O₁₂₈ kernel proved to be the best at t = 1, this is no longer true at t = 3 and N_smooth = 128, where the O₂₅₆ performs better across all resolutions at this time. And looking at the result at the earlier time of t = 1, the O₁₂₈ just barely performed better at N_smooth = 128 nx = 64, but worse than O₂₅₆ for higher resolution. We can also see that the O₂₅₆ performs better than the O₆₄ at N_smooth = 64 for the higher resolutions(n_x = 128, 256) at both t = 1 and t = 3. The improvement seen in N_smooth = 64 for GDSPH is the largest improvement that we see in across all our Gresho-Chan cases, where the optimized kernels effectively half the L₁,error compared to the Wendland C₂ and C₄ kernels, resulting in errors similar to N_smooth = 128 errors of the Wendland kernels (at both t = 1 and t = 3). This is quite clearly seen in Fig. 6, when comparing the O₆₄ N_smooth = 64 vortex structure with the N_smooth = 128, C₂ and C₄ vortex structure. A quite surprising result is that the C₄ kernel performs better than the C₂ kernel for N_smooth = 64, while performing worse for N_smooth = 32,128, where at N_smooth = 32 it performs drastically worse. For N_smooth = 32, we can see that the C₂ kernel performs relatively well, but with O₃₂ being about a 20-30% improvement at t = 1 and the only kernel with some sort of vortex structure at t = 3 (as can be seen in Fig. 6).

Looking at the result from ISPH, we can see a much less difference between all the different kernels, especially at high neighbor number (N_smooth = 128,256). At t = 3 and N_smooth = 256 we can see that the C₄ kernel actually performs better than all other kernels, including the optimized kernel O_I,256. This is quite surprising, as the O_I,256 performs better at earlier times, and looking at the E₀ errors at t = 3 in Fig. 9, we can see that they are higher for C₄ compared to all other kernels. This likely means that C₄ has smaller higher order errors, compared to the other kernels. The OI,₃₂ kernel performs significantly better than all the other kernels at N_smooth = 32 and is the only optimized kernel with negative Fourier transform at relatively low wave number, which means that it becomes unstable to pairing instability at higher neighbor numbers (above N_smooth = 64). The Wendland C₄ kernel is by far the worst kernel at lower neighbor numbers (N_smooth = 32,64) for ISPH. For N_smooth = 64 at both t = 1 and t = 3, the O_I,64 produce a significant improvement compared to the Wendland kernels.

The convergence scaling seems to be similar between kernels, but with some experiencing a flattening at higher resolution. The C₄ kernel seems to scale a bit better than other kernels at higher resolution, though still being worse in general. The scaling will also depend on the given dissipation scheme/parameters used. We have also run the Gresho-Chan vortex with different dissipation parameters (α = 0.01,0.05,0.1) and the relative results seen in this section remains the same.

Fig. 6 Azimuthal velocity in the radial direction of the Gresho-Chan vortex at t = 3. Here the resolution is nx = 64 and simulation was performed with GDSPH (left three columns) and ISPH (right three columns). We compare the optimized kernel at a given Nsmooth to the popular Wendland C2 and C4 kernels, and vary the number of neighbors in each row (Nsmooth = 32,64,128,256), referred here to N32, N64, N128, N256. The solid black line show the analytical solution and the black dots show the results from the simulation. The X axis (cylindrical radius) extends from 0.0 to 0.7, and the vertical extent from the top to the bottom of the analytical solution is 1. We can see a significant improvement in the optimized kernels for the GDSPH method for Nsmooth = 32,64,128, where the vortex structure becomes more pronounced. Overall less noise can be seen for all Nsmooth. The ISPH method captures the overall vortex structure well for all Nsmooth, due to exact linear gradients. But we can see significant improvement in the noise level of the optimized kernels, compared to the C2 and C4 kernels.

Fig. 7 Convergence study of the Gresho-Chan vortex with the GDSPH method for Nsmooth = 32,64,128,256. Here we plot the L1,error over resolution elements in the x direction (nx) at t = 3. Significant improvements can be seen for all optimized kernels compared to the Wendland C2 and C4 kernels. The kernels were optimized in respect to an earlier time t = 1 and we can see that O256 have an improved result over O128 for Nsmooth = 128 at this later time, even beating out O64 for Nsmooth = 64 at higher resolution (nx > 128). The C4 kernel exhibits slightly better convergence slope at high Nsmooth and nx > 128, likely due to it’s higher order properties (fewer O(h2) errors).

Fig. 8 Convergence study of the Gresho-Chan vortex with the ISPH method for Nsmooth = 32,64,128,256. Here we plot the L1,error over the resolution elements in the x direction (nx) at t = 3. Significant improvements can be seen for the optimized kernels at Nsmooth = 32,64 compared to the Wendland C2 and C4 kernels. While at higher Nsmooth, there is much less difference between all the kernels. The kernels were optimized in respect to an earlier time, t = 1, and we observed that C4 performs better at this later time for Nsmooth = 256 than all the optimized kernels at all resolutions. This is surprising, as E0 errors of the optimized kernels are much lower than the C4 kernel (see Fig. 9), which might indicate that long-term behavior at this Nsmooth is heavily determined by the higher order errors (O(h2)).

Fig. 9 Convergence study of the Gresho-Chan vortex with the ISPH method for Nsmooth = 256. Here we plot the E0 error versus the number of resolution elements in the x direction (nx) at t = 3. Much lower zeroth order errors for the optimized kernel and C2 kernel can be seen compared to the higher order C4 kernel. Note here that even though being zeroth order (dependent on local particle order), the error still decreases with increasing resolution ($n_x^{-0.55}$) because the particle distribution becomes more uniform at higher resolution.

4.3 Kelvin Helmholtz

The Kelvin-Helmholtz (KH) instability occurs when there is a shearing flow between two fluid interfaces. The instability quickly generates a series of rolling waves at the interface between the two fluids, which eventually breaks down into turbulence, mixing the two fluids. The instability plays a significant role in all kinds of fluids and is a crucial component for any hydrodynamical code to model, as it effectively allows two fluids to mix and generate turbulence. However, traditional SPH methods have been known to struggle with triggering the KH instability, particularly when there are significant density contrasts between the fluids. In this paper we model the KH instability using a similar setup to the ones outlined in McNally et al. (2012); Lecoanet et al. (2016); Tricco (2019), using a thin periodic box ($L=[1,1,2\frac{\sqrt{6}}{n_x}]$), here n_x is the average number of particles in the x direction, and the z boundary is set to be roughly 24 particle spacings. A smoothing function was applied to the density and the shearing velocity accordingly: $$f_R=0.5\ast (\tanh \frac{y+y_{inner}}{\delta }-\tanh \frac{y-y_{inner}}{\delta }),$$(51) $$\rho ={\rho }_{outer}+f_R({\rho }_{inner}-{\rho }_{outer}),$$(52) $$v_x=v_{x,outer}+f_R(v_{x,inner}-v_{x,outer}),$$(53) $$v_y=0.01\sin (4\pi x).$$(54)

Here we set ρ_outer = 1, δ = 0.025, v_x,outer = −0.5, and v_x,inner = 0.5. We simulated the KH instability with two different density contrasts (ρ_inner = ∆ρ = 2 and ρ_inner = ∆ρ = 8). The pressure was initially uniform ($P_0=\frac{1}{\mathrm{\Gamma }{\mathcal{M}}^2}$), where Γ is the adiabatic index and M is the Mach number. We used an adiabatic equation of state with adiabatic index of Γ = 5/3 and a Mach number of $\mathcal{M}=\sqrt{6}/5$. The IC is setup using the IC glass generator as described in Sect. 3.2 for all the different kernels. We ran the simulation until t = 4 for each of the optimized kernels at their specific N_smooth = 32, 64, 128, 256 and SPH method. We then compared these kernels with the Wendland C₂ and C₄ kernel.

In Fig. 10 we show the rendered structure of the KH instability for the GDSPH ∆ρ = 2 and ∆ρ = 8 at t = 2 and t = 4. Due to the N_smooth = 256 simulations producing near identical results to the N_smooth = 128 we only show the results of N_smooth = 128. For both the low and high density contrast cases, we can see a clear improvement in the growth of the instability for the optimized kernels, producing similar structure as the C₄ kernel at half the N_smooth. In Fig. 12 we show the amplitude growth of the k = 2π mode, using the same procedure as outlined in McNally et al. (2012), which more quantitatively show the improvement of the optimized kernel for GDSPH.

In general, the difference between the kernels in the ISPH method is very minor for the KH case, across the different N_smooth, density contrasts and times. The only exception is for N_smooth = 32 at t = 4 (see Fig. 11), where in the ∆ρ = 2 case the C₄ experiences significant non-conservation of angular momentum, leading to the flow rotating (still capturing the rolling vortex waves). For ∆ρ = 8, t = 4 the optimized kernel, produces a vortex structure (seen at higher N_smooth for all kernels), while the C₄ kernel fails to capture the vortex. The optimized kernel also produces sharper density structure compared to C₂ kernel.

We have also performed more comparisons and a resolution study between GDSPH and ISPH for the KH test, which can be found in Appendix A.

Fig. 10 Density rendering of the Kelvin-Helmholtz instability at t=2 (left column) and t=4 (right column) for an initial density contrast of ∆ρ = 2 (top panels) and ∆ρ = 8 (bottom panels) with the GDSPH method. The X and Y axes have a length of one. Here one can see a comparison in the development of the instability between the optimized kernel at a given Nsmooth and the popular Wendland C2 and C4 kernels. We observed a clear improvement in the optimized kernel, reproducing a similar behavior of the C4 kernel at double the neighbor count (Nsmooth). For the high density contrast ∆ρ = 8, there are only minor differences between the different kernels at Nsmooth = 128. These can be compared to the higher resolution cases given in Appendix A.

Fig. 11 Density rendering of the Kelvin-Helmholtz instability at t = 4 for Nsmooth = 32 with the ISPH method. The X and Y axes have a length of one. Here one can see a comparison in the development of the instability between the optimized kernel (OI, 32) and the popular Wendland C2 and C4 kernels. We have selected three cases where discernible differences can be seen (nx = 128, ∆ρ = 2; nx = 64, ∆ρ = 2; nx = 128, ∆ρ = 8). For the ∆ρ = 2 cases we observed that the C2 and C4 kernels exhibit more angular momentum errors (inherent to ISPH) compared to the OI,32 kernel, where the inner stream becomes significantly twisted, deviating from the high resolution and high Nsmooth results (see Appendix A). For the high density case nx = 128, ∆ρ = 8, we observed that C4 is unable to generate the same vortex structure in the upper plane of the stream and exhibit some asymmetry around the plane.

4.4 Sod shock tube

Here we test the classic Sod shock tube (Sod 1978), where a fluid medium is initially divided into two separate regions with different pressures and densities. The initial conditions are ρ_left = 1, P_left = 1 for x < 0 and ρ_right = 0.25, P_right = 0.1795 for x > 0. We use an adiabatic equation of state with adiabatic index of Γ = 5/3. At the start of the simulation, we imagine removing a wall between the two regions creating a discontinuous state between the fluid’s generating a shockwave. The produced shockwave has three distinct regions, the right-moving shock front, the left-moving rarefaction wave, and the contact discontinuity. Here there are once again many ways to setup the initial distribution of particles. This case is particularly difficult due to the discontinuous density profile, as it is not possible to follow a discontinuous solution with the smooth density estimate. That is, even if two different density glasses or lattices are set up, the density should still be smoothed at the interface. One could either try one’s best to relax a glass following the discontinuous solution or smooth the initial shock by a tiny amount (proportional to the smoothing length). An argument against the complete discontinuous solution (in density) would be that it would rarely happen in a dynamic SPH simulation, as it would always be accompanied with some smoothing. Other quantities, such as velocities and thermal energy, of course occur discontinuously in simulations. We chose to smooth the discontinuity on the scale of the smoothing length. We used a thin periodic 3D box to approximate the 1D test, with L_x = 2.0 and L_y = L_z set to be roughly 16 particle spacings in the low density region. We used L_x = 2.0 due to using periodic boundary conditions such that the shock generated at the periodic boundary does not interact with the central shock. We used a resolution of around n_x = 64 × 2 in the initial high density region (the factor of two is to account for the box extension). We used a constant α = 1 and β = 2 for the artificial viscosity.

As the velocity shows the biggest difference, we highlight the velocity structure between the different kernels in Fig. 13 for N_smooth = 64. From this figure we can see that the optimized kernel’s exhibit less noise in the post-shock region compared to the Wendland C2 and C4 kernels for both GDSPH and ISPH. Other than that, there is not much difference between the kernels in this case.

Fig. 12 Amplitude of the k = 2π mode over time for different smoothing kernels and density contrasts: ∆ρ = 2 (top) and ∆ρ = 8 (bottom). The ISPH method shows little difference between the different smoothing kernels, capturing a growth rate of roughly γ ∝ π in ∆ρ = 2 case and a slightly lower growth rate for ∆ρ = 8 (as expected). The GDSPH method of Nsmooth = 64 struggles to capture the same growth rate as the ISPH method, but we observed that the optimized kernel improves the growth rate of the KH mode for both ∆ρ = 2 and ∆ρ = 8.

Fig. 13 Velocity structure of the Sod-shocktube test at t = 0.2 performed at nx = 64 and Nsmooth = 64. The thin black line shows the exact solution. One can see that the optimized kernel exhibit less noise in the post-shock region compared to the C2 and C4 kernel.

4.5 Evrard collapse

The Evrard collapse test (Evrard 1988; Steinmetz & Mueller 1993; Wadsley et al. 2017; Mandal et al. 2023), simulates the collapse of adiabatic gas, which include a high Mach shock and a cold pre-shock infall. The density is defined as $$\rho (r)=\frac{M(R)}{2\pi R^{2}r},$$(55)

where the initial cloud radius is R = 1 and the initial cloud mass is M = 1. The internal energy is u = 0.05GM/R and we use an adiabatic index of γ = 5/3. We setup the initial density using our IC generator, using roughly 28 000 particles. As this is an open boundary IC, we need to assume an initial low density medium for the IC generator to fit against. This is because the IC generator uses only the SPH forces to fit the density profile, this requires a periodic boundary and density to be defined everywhere. After the IC generator has fitted the density profile, we remove all particles beyond R = 1 and ensure that particle mass sum up to M = 1.

The biggest difference in the Evrard collapse (at t = 0.08) between the two methods GDSPH and ISPH and the different kernels can be seen in the noise level, similar to the Sod shocktube test. In Fig. 14 we can see the velocity structure of the two SPH methods and kernels. The noise level of C₄ is slightly larger than the C₂ and the optimized kernels. With C₄ having a slightly sharper shock front than the other two.

4.6 Subsonic blob test

The blob test consist out of a high density spherical cloud at rest within a uniformly flowing low density medium. As the high density blob is accelerated by the surrounding medium, fluid instabilities such as Kelvin-Helmholtz and Rayleigh-Taylor instabilities disrupts the cloud and mixes it with the surrounding medium. The original blob test of Agertz et al. (2007), models a wind with a super sonic Mach number of M = 2.7. However, we have seen from the Sod shock tube and Evrard collapse tests, that the quality of the smoothing kernel does not have much effect on high Mach tests. Particle noise and kernel errors are more relevant for subsonic flows, and as such we decided to instead model a subsonic blob test with a Mach number of M = 0.5 to test the effect of the kernels on the mixing in the subsonic regime. The domain is periodic, with dimensions (10,10,40) in units of cloud radius. The cloud was put in pressure equilibrium and has a density ratio of Pρ = ρ_cloud/ρ_wind = 100. The wind velocity was set by the Mach number v_wind = Mc_s, and we defined the characteristic cloud crushing time to be $t_{crush}=\sqrt{\mathrm{\Delta }\rho }{R}_{cloud}/v_{wind}$ and ran the simulation for 10t_crush. In our analysis of the blob test, we defined the cloud mass to be M_cloud = M(ρ > ρ_cloud/3) and the intermediate temperature mass as in Braspenning et al. (2023); Sandnes et al. (2025): $$\begin{array}{r}{M}_{mix}=\sum m_i(\log (T_{mix})-0.25\log (\mathrm{\Delta }\rho )<\log T_i\\ <\log (T_{mix})+0.25\log (\mathrm{\Delta }\rho )).\end{array}$$(56)

Here, m_i is mass of the particle, T_i is the temperature of the particle, and T_mix is the geometric mean temperature: $T_{mix}=\sqrt{T_{cloud}{T}_{wind}}=T_{cloud}\sqrt{\mathrm{\Delta }\rho }$.

From the top panel of Fig. 16, we can see that all models disrupts the cloud and reduces the high density mass in a few cloud crushing times (3-5 t_crush). However, it is clear from the rendering of the simulation at t = 5t_crush in Fig. 15 that the mixing of the gas with the surrounding medium differs between the different SPH methods and smoothing kernels. This can also be seen in the intermediate temperature mass evolution, where initially this mass increase is similar in all models, but the subsequent mixing of this mass with the warm medium evolves very differently. ISPH mixes the cloud more aggressively than GDSPH, where for ISPH at t = 5t_crush most of the intermediate temperature mass have mixed with the surrounding medium, while GDSPH have majority of the initial cloud mass in this intermediate temperature. The result of the intermediate temperature mass, is reminiscent of the KH mode amplitude growth dependency that we saw in Fig. 12.

Fig. 14 Velocity structure of the Evrard collapse at t = 0.8 performed with 28000 particles and Nsmooth = 64. The thin black line shows the reference solution from a 1D PPM simulation with 8192 grid points (Mandal et al. 2023). One can see that the optimized kernel and the C2 kernel exhibit less noise in the post-shock region compared to the C4 kernel.

Fig. 15 Density rendering of the subsonic blob test (∆ρ = 100, M = 0.5 Rcloud = 0.1) at t = 5tcrush with a resolution of nx = 64 in the low density medium. The three left columns depict the cloud state with the ISPH method and the kernels (OI,64, C2, C4), and three right columns show the GDSPH method and the kernels (O64, C2, C4). The X axis has a length of one, and Y axis has a length of three. One can see that the ISPH method disrupts and mixes the cloud more aggressively than the GDSPH method. We also observed that we increase mixing with the optimized kernels compared to the C2 and C4 kernel in both the GDSPH and ISPH method.

Fig. 16 Time evolution of the cloud mass ($M(\rho >\frac{1}{3}{\rho }_{cloud})$; top column) and the mass of intermediate temperature gas (Eq. (56); bottom column). Once can see that it takes about the same time for the cloud mass to be lost in all the different modules, though with a slight faster decrease in the ISPH modules, with the fastest decrease being with the OI,64 kernel + ISPH. However, the amount of mixing in intermediate temperature gas differs quite a bit, where the optimized kernels increases the amount of mixing in both the ISPH and GDSPH models. The mixing results are reminiscent of the KH mode growth results shown in Fig. 12.

5 Discussion

In this paper we have generated a set of new optimized kernels through the linear combination of other positive-definite kernels whose coefficients were optimized for the GDSPH method for the hydrostatic glass, the Gresho-Chan vortex, and for the ISPH method for the Gresho-Chan vortex. We find that the kernels optimized for the hydrostatic glass perform significantly worse for dynamical cases. We find that all kernels that had regular large oscillations in their Fourier transform performed worse for the Gresho-Chan vortex. The kernels produced by optimizing to the Gresho-Chan vortex all produce a quite steady decreasing Fourier transform with very slight oscillatory behavior and with a very similar functional form to that of the C₂ kernel. However, compared to the C₂ kernel, the majority of optimized kernels have a Fourier transform that decreases more rapidly with wave numbers (past k > 3π). We observed that in general, the Fourier transform decreases faster past k > 3π for the kernels optimized at a higher neighbor number, and it decreases faster before k < 3π for the kernels optimized at a lower neighbor number. Having a lower power at a high wave number seems to be important for higher neighbor numbers, while a lower power at low wave numbers is important for smaller neighbor numbers.

The new optimized kernels produce significant improvements for the GDSPH method, showing an improvement for the Gresho-Chan vortex (n_x = 64 t = 1) L_1,error of −20, −41, −31, and −17% relative to the L_1,error of the C₂ kernel for N_smooth = 32,64,128,256, respectively. The optimized kernel for N_smooth = 64 gives roughly the same L_1,error as that of N_smooth = 128 for the C₂ kernel. At the later time, we observed that the O₂₅₆ kernel performs better than O₁₂₈ for all n_x and even better than O₆₄ at higher resolution (n_x = 128,256), as O₂₅₆ performs quite well even at t = 1 for both N_smooth = 64 and N_smooth = 128. One could use O₂₅₆ for all N_smooth > 64. The C₄ kernel has the steepest improvement for N_smooth = 256, remaining near linear with resolution, while O₂₅₆ starts to flatten out beyond n_x = 256. This is likely due to the fact that O₂₅₆ reduces zero- and first-order errors a lot, while C₄ reduces higher-order errors more than O₂₅₆. For ISPH, we also observed improvements of the L_1,error in the Greso-Chan vortex, which are −30, −21, −11, and −8% relative to the L_1,error of the C₂ kernel for N_smooth = 32,64,128,256, respectively. This is visible mainly in a reduction of noise. At the later time, there is much less difference between the kernels at the higher N_smooth = 128,256, while a significant improvement can still be seen for N_smooth = 32, 64. The C₂ kernel is in general the best-performing smoothing kernel from the non-optimized kernels for most N_smooth. At later times, at t = 3, we observed that C₄ performs better at N_smooth = 64 for GDSPH, the reason for which is not fully clear, and N_smooth = 256 for the ISPH method, where C₄ is the best-performing kernel of all at t = 3, likely due to the leading orders being much lower and thus giving more importance to higher-order errors. The noise of N_smooth = 256 for the ISPH method for C₄ is still greater than the O_I,256 kernel, as can be seen in Fig. 6. But the peak of the vortex is slightly higher with the C₄ kernel. We can compare the convergence rate ($n_x^{-p}$) to Cabezón & García-Senz (2024), who used ISPH and mix Sinc kernels to provide a more accurate and stable kernel against pairing instability. Those authors show a convergence rate of p = 0.2 at N_smooth = 60 to p = 1.03 at N_smooth = 400 after t = 1. For ISPH/GDSPH(t = 3)[t = 1] and the optimized kernels, we show a convergence rate of (p = 0.5/p = 0.3)[p = 0.4/p = 0.3] at N_smooth = 32, /(p = 0.6/p = 0.5)[p = 0.5/p = 0.5] at N_smooth = 64, (p = 0.9/p = 0.7)[p = 0.8/p = 0.7] at N_smooth = 128 to (p = 1.1/p = 0.9)[p = 1.0/p = 0.9] at N_smooth = 256 after (t = 3)[t = 1]. In Wadsley et al. (2017), with modern artificial viscosity switches, they found a convergence rate of p = 0.8 for N_smooth = 200 with the C₄ kernel at t = 3. The moving mesh code AREPO in Springel (2010) yields a convergence rate of p = 1.4. Note that comparison to other modern SPH codes is slightly unfair, as we used a classic artificial viscosity scheme with constant α and β values, where we have chosen a low α value that gives a good result for the Gresho-Chan vortex.

For the Kelvin-Helmholtz instability, the biggest improvement could be seen for the GDSPH method, where the new optimized kernels captured similar results to that of the C₄ kernel at double N_smooth (both at ∆ρ = 2 and ∆ρ = 8 density contrast). For the ISPH method, the main improvement could be seen in angular momentum conservation of the low N_smooth runs, where O_I,32 performs better than C₂ and C₄. For the Sod-shocktube and Evrard collapse test, the main improvement of the optimized kernels could be seen as a reduction of the noise in the post-shock region.

In conclusion, we find that the linear-combination of kernels generates optimized kernels that significantly improve results, particularly in the low N_smooth regime. This can give a result similar to doubling the N_smooth of previous smoothing kernels (C₂, C₄) without any additional computational cost.

For ISPH, the results of the linear fluid behavior is only weakly dependent on the smoothing kernel and neighbor number (as seen in the Gresho-Chan vortex and the KH instability). We could see that the particle noise, angular momentum errors, and the mixing of the subsonic blob test improved with the new optimized kernels. For GDSPH, we observed an overall improvement with the optimized kernels for all test cases. From these tests we observed that ISPH give a significant improvement in convergence and behavior for subsonic flows compared to GDSPH. Results from the high Mach tests such as the Sod shocktube and Evrard collapse are much more insensitive to the SPH method and smoothing kernel used. For GDSPH, we recommend using O₃₂ for N_smooth <= 32, O₆₄ for N_smooth <= 96, and O₂₅₆ for higher values of N_smooth. For ISPH, we recommend using O_I,32 for N_smooth <= 32, O_I,64 for N_smooth <= 96, and O_I,256 for higher values of N_smooth.

We have introduced a new IC generator to generate glass distributions of arbitrary density profiles in SPH. We have also introduced self-bias corrections to the kernel, that have been adjusted to the glass distribution instead of the lattice distribution of earlier work. This ensures more accurate density interpolation for the glass distribution of positive-definite kernels used in this paper.

This work shows that targeted kernel optimization can enhance SPH convergence beyond what is achievable with the standard compact kernels commonly used. It also demonstrates that the most effective kernel for a given problem depends on the specific SPH formulation used. These findings motivate future exploration of new optimization strategies and further improvements to both the smoothing kernels and the SPH method as a whole.

Data availability

The initial condition generator, optimized kernels, and test cases used in this work are publicly available at: https://github.com/robertwissing/testsuite.

Acknowledgements

We acknowledge support from the European High Performance Comput- ing Joint Undertaking (EuroHPC JU) and the Research Council of Norway through the funding of the SPACE Center of Excellence (grant agreement N0 101093441). Parts of the simulations were performed using the resources from the National Infrastructure for High Performance Computing and Data Storage in Norway, UNINETT Sigma2, allocated to Project NN9477K. We acknowledge the use of SPLASH (Price 2007) for the visualization of the results.

References

Appendix A Differences between GDSPH and ISPH

In this appendix we include a resolution study between GDSPH and ISPH for the KH instability using the new optimized kernels for each respective N_smooth, which can be seen in Fig. A.1 for ∆ρ = 2 and Fig. A.2 for ∆ρ = 8. From these figures it is clear that the result for ISPH is highly independent of N_smooth, producing the same solution for N_smooth = 32 and N_smooth = 128. The GDSPH results are on the other hand, much more dependent on N_smooth, and exhibit severe damping of the instability at N_smooth = 32. At N_smooth = 128, the GDSPH method produces similar structures to that of the ISPH method at ∆ρ = 2. For high density contrast ∆ρ = 8, we can see that high N_smooth for the ISPH method generates more noise at the boundary layer at t = 2 compared to GDSPH in this case. We can also see that at later time, the high density structures are more dissipate in the ISPH method compared to the GDSPH method. We still see much more consistent behavior in the ISPH method compared to GDSPH method.

The high density contrast KH instability, is interesting as it is usually accompanied with strong surface tension effects in the Lagrangian formalism of SPH, due to high E₀ and E₁ errors at the discontinuous boundary. In GDSPH, these errors are reduced by doing a geometric density weighting in the gradient operator, reducing the errors at the density contrast. This sort of correction, deviates from the density estimate derivation of Lagrangian SPH and will thus introduce an entropy error. In GDSPH this is corrected for so that the method remains second-order accurate in entropy. ISPH removes the linear gradient error from it’s gradients, which effectively removes more of the "partition" force that exist in SPH and only have E₀ error forces to remain ordered. This sort of correction also introduces entropy errors, which would potentially be first order in this case. In Fig. A.3, we show how one of the KH swirls develops when not using diffusion versus using diffusion for both the GDSPH and ISPH methods. Significant amounts of entropy bubbles can be seen in the ISPH version without diffusion, and a few can be seen in GDSPH. Adding thermal diffusion (Shen et al. 2010) allows for local mixing between the cold and hot phases, and we can see that we get good behavior in both ISPH and GDSPH. Slight numerical surface tension effects can be seen for GDSPH, while none can be seen for ISPH.

Appendix B Kernel self-bias correction coefficients

In this section we present the self-correction bias coefficients present in $${\rho }_a={\rho }_{a,estimate}-ϵm_{a}W(0,h_a)$$(B.1) $$ϵ(N_{smooth})=\frac{c_{i+1}-c_i}{N_{min,i+1}-N_{min,i}}(N_{smooth}-N_{min,i})+c_i,$$(B.2)

where N_min,i ∈ {16,24,32,48,64,96,128,256,512). The c_i coefficients can be found in Table B.1 for all kernels.

Appendix C Kernel tables

The linear-combination of kernels after optimization are fitted are fitted to a polynomial in the form $$P(x)=\sum _{i=0}^8{c}_i{x}^i,$$(C.1) $$W(q)=f_{norm}(1+q^{2}P(q)).$$(C.2)

The polynomial coefficients for all the optimized kernels are given in C.1 and C.2.

Appendix D Code implementation

A simple code snippet of the kernel function is included for clarification and for easier implementation, for both GDSPH and ISPH.

Fig. A.1 Density rendering of the Kelvin-Helmholtz instability (∆ρ = 2) at t = 2 (three left columns) and t = 4 (three right columns) for varying resolution (nx = 64,128,256), Nsmooth = 32,64,128, and method (GDSPH upper three rows and ISPH bottom three rows). The X and Y axis have a length of 1a. These simulations were all run with the optimized kernel for each ONsmooth GDSPH and OI,Nsmooth for ISPH. We can see that ISPH manages to generate the same vortex structure, regardless of the number of neighbors if the resolution is high enough, while GDSPH struggles with lower number of neighbors even at the highest resolution.

Fig. A.2 Density rendering of the Kelvin-Helmholtz instability (∆ρ = 8) at t = 2 (three left columns) and t = 4 (three right columns) for varying resolution (nx = 64,128,256), Nsmooth = 32,64,128, and method (GDSPH upper three rows and ISPH bottom three rows). The X and Y axis have a length of 1. These simulations were all run with the optimized kernel for each ONsmooth GDSPH and OI,Nsmooth for ISPH. Similar as the low density version of Fig. A.1, we can see that ISPH captures instability well regardless of neighbor number. We can see additional diffusion from the ISPH runs compared to the GDSPH runs, even with the highest resolution. There also seem to be additional noise in the small-scale density structures as we increase the number of neighbors for ISPH. Potentially related to the issue found in Fig. A.3

Fig. A.3 Particle density rendering of Kelvin-Helmholtz instability swirl at t = 2 for ∆ρ = 8. The figures shows the different behavior of not using diffusion versus using diffusion for the GDSPH and ISPH methods. A significant number of entropy clumps can be seen being generated in the low density medium for ISPH, as the removal of linear gradient together with sharp boundaries causes chaotic noise at the boundary. Lagrangian SPH formalisms that are derived directly from the density estimate, provide full spatial conservation of entropy. However, deviation from the Lagrangian formalism, while still using the traditional density estimate means that we will introduce entropy errors in our solution. GDSPH corrects for firstorder errors in entropy making it second-order in entropy. ISPH removes linear-gradient errors and will thus only be first-order accurate with entropy. The issue might be further aggravated by the use of average gradient kernels within ISPH. Adding thermal diffusion allows for local mixing between the cold and hot phase, and we can see that we get good behavior in both ISPH and GDSPH. Slight numerical surface tension effect can be seen for GDSPH, while none can be seen for ISPH.

All Tables

All Figures

Fig. 1 Relaxed glass IC generated by the IC generator outlined in Sect. 3.2. Here we model the IC of the Kelvin Helmholtz instability (Sect. 4.3), with a density contrast of ∆ρ = 8. Top : particle distribution within a thin slice of thickness ∆ z = 0.015. Bottom : density profile in the y direction. The red line depicts the analytical density profile and the black dots the density of the particles. The glass generated density profile is in excellent agreement with the analytical solution.

In the text

Fig. 2 Top : optimized kernel functions for a glass distribution (Sect. 3.4.1) together with the missing Wendland kernels (CM04, CM05), the Buhmann kernel (BUH), and the Wu kernel (WU2). Middle : corresponding derivative of these kernel function. Bottom: fourier transform of these kernels, scaled by k/10. The optimized kernels are quite similar, with some small variations being seen in the Fourier transform. The Oglass,32 kernel deviates the most from the other kernel, with less deep troughs in the Fourier transform. These kernel’s generate excellent glass distributions but degrade significantly in dynamical simulations.

In the text

Fig. 3 Top : optimized kernel functions for the Gresho-Chan vortex (Sect. 3.4.2), for both GDSPH and ISPH together with the popular Wendland kernels (C2, C4). Middle : corresponding derivative of these kernel function. Bottom : fourier transform of these kernels, scaled by k/10. The optimized kernel exhibit a clear trend in it’s Fourier transform, where the kernels optimized for higher Nsmooth emphasizes lower power beyond k > 3π. The general form of the optimized kernel is quite similar to the Wendland C2 kernel. In general GDSPH produces lower power than ISPH in the high wave number regime for a given Nsmooth. The OI,32 kernel is the only optimized kernel that exhibits a negative Fourier transform at low wavenumbers. Among the other optimized kernels, only O256 has a negative Fourier transform, but this occurs at higher wavenumbers and remains stable to pairing instability up to at least Nsmooth < 500.

In the text

Fig. 4 Relative increase or decrease of the δρerr (top), E0 (middle), and E1 (bottom) of the hydrostatic glass when relaxed for different kernels in respect to the corresponding error of the C2 kernel. The simulations here were performed with only the GDSPH method. The terms N32, N64, N128, and N256 refer to the number of neighbors: Nsmooth = 32,64,128,256, respectively. The δρerr, E0 and E1 of the C2 kernel for each Nsmooth is δρerr = 5.539 × 10−3,3.747 × 10−3,1.243 × 10−3,5.860 × 10−4, E0 = 2.763 × 10−3,2.024 × 10−3, 6.797 × 10−4,4.253 × 10−4, and E1 = 2.210 × 10−2,1.484 × 10−2,2.557 × 10−3,8.305 × 10−4. The optimized kernels can in general be seen to drastically reduce the E1 error relative to the C2 kernel. For the glass optimized kernel’s we can also see a general decrease in both δρ and E0 for their given Nsmooth.

In the text

Fig. 5 Relative increase or decrease of the L1,error of the Gresho-chan vortex at t = 1 and nx = 64 for different kernels in respect to the L1,error of the C2 kernel. The terms N32, N64, N128, and N256 refer to the number of neighbors: Nsmooth = 32,64,128,256. Left : simulations done with the GDSPH method. Right : simulations done with the ISPH method. The L1,error of the C2 kernel for each Nsmooth is L1,error = 1.49375 × 10−1, 9.57342 × 10−2, 4.70888 × 10−2, 2.99816 × 10−2 for GDSPH and L1,error = 9.34635 × 10−2, 4.66192 × 10−2, 2.64744 × 10−2, 2.22268 × 10−2 for ISPH. An improvement for the optimized kernel at a given Nsmooth across all Nsmooth can be seen for both methods, where the biggest improvements are for the N64 and N128 for GDSPH and N32 and N64 for ISPH. Of the non-optimized kernels, we observed that the C2 kernel performs best, except for Nsmooth = 32 (GDSPH), Nsmooth = 32,64 (ISPH), where CM05 performs better and C4 performs slightly better at Nsmooth = 256 (ISPH).

In the text

Fig. 6 Azimuthal velocity in the radial direction of the Gresho-Chan vortex at t = 3. Here the resolution is nx = 64 and simulation was performed with GDSPH (left three columns) and ISPH (right three columns). We compare the optimized kernel at a given Nsmooth to the popular Wendland C2 and C4 kernels, and vary the number of neighbors in each row (Nsmooth = 32,64,128,256), referred here to N32, N64, N128, N256. The solid black line show the analytical solution and the black dots show the results from the simulation. The X axis (cylindrical radius) extends from 0.0 to 0.7, and the vertical extent from the top to the bottom of the analytical solution is 1. We can see a significant improvement in the optimized kernels for the GDSPH method for Nsmooth = 32,64,128, where the vortex structure becomes more pronounced. Overall less noise can be seen for all Nsmooth. The ISPH method captures the overall vortex structure well for all Nsmooth, due to exact linear gradients. But we can see significant improvement in the noise level of the optimized kernels, compared to the C2 and C4 kernels.

In the text

Fig. 7 Convergence study of the Gresho-Chan vortex with the GDSPH method for Nsmooth = 32,64,128,256. Here we plot the L1,error over resolution elements in the x direction (nx) at t = 3. Significant improvements can be seen for all optimized kernels compared to the Wendland C2 and C4 kernels. The kernels were optimized in respect to an earlier time t = 1 and we can see that O256 have an improved result over O128 for Nsmooth = 128 at this later time, even beating out O64 for Nsmooth = 64 at higher resolution (nx > 128). The C4 kernel exhibits slightly better convergence slope at high Nsmooth and nx > 128, likely due to it’s higher order properties (fewer O(h2) errors).

In the text

Fig. 8 Convergence study of the Gresho-Chan vortex with the ISPH method for Nsmooth = 32,64,128,256. Here we plot the L1,error over the resolution elements in the x direction (nx) at t = 3. Significant improvements can be seen for the optimized kernels at Nsmooth = 32,64 compared to the Wendland C2 and C4 kernels. While at higher Nsmooth, there is much less difference between all the kernels. The kernels were optimized in respect to an earlier time, t = 1, and we observed that C4 performs better at this later time for Nsmooth = 256 than all the optimized kernels at all resolutions. This is surprising, as E0 errors of the optimized kernels are much lower than the C4 kernel (see Fig. 9), which might indicate that long-term behavior at this Nsmooth is heavily determined by the higher order errors (O(h2)).

In the text

Fig. 9 Convergence study of the Gresho-Chan vortex with the ISPH method for Nsmooth = 256. Here we plot the E0 error versus the number of resolution elements in the x direction (nx) at t = 3. Much lower zeroth order errors for the optimized kernel and C2 kernel can be seen compared to the higher order C4 kernel. Note here that even though being zeroth order (dependent on local particle order), the error still decreases with increasing resolution ($n_x^{-0.55}$) because the particle distribution becomes more uniform at higher resolution.

In the text

Fig. 10 Density rendering of the Kelvin-Helmholtz instability at t=2 (left column) and t=4 (right column) for an initial density contrast of ∆ρ = 2 (top panels) and ∆ρ = 8 (bottom panels) with the GDSPH method. The X and Y axes have a length of one. Here one can see a comparison in the development of the instability between the optimized kernel at a given Nsmooth and the popular Wendland C2 and C4 kernels. We observed a clear improvement in the optimized kernel, reproducing a similar behavior of the C4 kernel at double the neighbor count (Nsmooth). For the high density contrast ∆ρ = 8, there are only minor differences between the different kernels at Nsmooth = 128. These can be compared to the higher resolution cases given in Appendix A.

In the text

Fig. 11 Density rendering of the Kelvin-Helmholtz instability at t = 4 for Nsmooth = 32 with the ISPH method. The X and Y axes have a length of one. Here one can see a comparison in the development of the instability between the optimized kernel (OI, 32) and the popular Wendland C2 and C4 kernels. We have selected three cases where discernible differences can be seen (nx = 128, ∆ρ = 2; nx = 64, ∆ρ = 2; nx = 128, ∆ρ = 8). For the ∆ρ = 2 cases we observed that the C2 and C4 kernels exhibit more angular momentum errors (inherent to ISPH) compared to the OI,32 kernel, where the inner stream becomes significantly twisted, deviating from the high resolution and high Nsmooth results (see Appendix A). For the high density case nx = 128, ∆ρ = 8, we observed that C4 is unable to generate the same vortex structure in the upper plane of the stream and exhibit some asymmetry around the plane.

In the text

Fig. 12 Amplitude of the k = 2π mode over time for different smoothing kernels and density contrasts: ∆ρ = 2 (top) and ∆ρ = 8 (bottom). The ISPH method shows little difference between the different smoothing kernels, capturing a growth rate of roughly γ ∝ π in ∆ρ = 2 case and a slightly lower growth rate for ∆ρ = 8 (as expected). The GDSPH method of Nsmooth = 64 struggles to capture the same growth rate as the ISPH method, but we observed that the optimized kernel improves the growth rate of the KH mode for both ∆ρ = 2 and ∆ρ = 8.

In the text

	Fig. 13 Velocity structure of the Sod-shocktube test at t = 0.2 performed at nx = 64 and Nsmooth = 64. The thin black line shows the exact solution. One can see that the optimized kernel exhibit less noise in the post-shock region compared to the C2 and C4 kernel.
In the text

	Fig. 14 Velocity structure of the Evrard collapse at t = 0.8 performed with 28000 particles and Nsmooth = 64. The thin black line shows the reference solution from a 1D PPM simulation with 8192 grid points (Mandal et al. 2023). One can see that the optimized kernel and the C2 kernel exhibit less noise in the post-shock region compared to the C4 kernel.
In the text

Fig. 15 Density rendering of the subsonic blob test (∆ρ = 100, M = 0.5 Rcloud = 0.1) at t = 5tcrush with a resolution of nx = 64 in the low density medium. The three left columns depict the cloud state with the ISPH method and the kernels (OI,64, C2, C4), and three right columns show the GDSPH method and the kernels (O64, C2, C4). The X axis has a length of one, and Y axis has a length of three. One can see that the ISPH method disrupts and mixes the cloud more aggressively than the GDSPH method. We also observed that we increase mixing with the optimized kernels compared to the C2 and C4 kernel in both the GDSPH and ISPH method.

In the text

Fig. 16 Time evolution of the cloud mass ($M(\rho >\frac{1}{3}{\rho }_{cloud})$; top column) and the mass of intermediate temperature gas (Eq. (56); bottom column). Once can see that it takes about the same time for the cloud mass to be lost in all the different modules, though with a slight faster decrease in the ISPH modules, with the fastest decrease being with the OI,64 kernel + ISPH. However, the amount of mixing in intermediate temperature gas differs quite a bit, where the optimized kernels increases the amount of mixing in both the ISPH and GDSPH models. The mixing results are reminiscent of the KH mode growth results shown in Fig. 12.

In the text

Fig. A.1 Density rendering of the Kelvin-Helmholtz instability (∆ρ = 2) at t = 2 (three left columns) and t = 4 (three right columns) for varying resolution (nx = 64,128,256), Nsmooth = 32,64,128, and method (GDSPH upper three rows and ISPH bottom three rows). The X and Y axis have a length of 1a. These simulations were all run with the optimized kernel for each ONsmooth GDSPH and OI,Nsmooth for ISPH. We can see that ISPH manages to generate the same vortex structure, regardless of the number of neighbors if the resolution is high enough, while GDSPH struggles with lower number of neighbors even at the highest resolution.

In the text

Fig. A.2 Density rendering of the Kelvin-Helmholtz instability (∆ρ = 8) at t = 2 (three left columns) and t = 4 (three right columns) for varying resolution (nx = 64,128,256), Nsmooth = 32,64,128, and method (GDSPH upper three rows and ISPH bottom three rows). The X and Y axis have a length of 1. These simulations were all run with the optimized kernel for each ONsmooth GDSPH and OI,Nsmooth for ISPH. Similar as the low density version of Fig. A.1, we can see that ISPH captures instability well regardless of neighbor number. We can see additional diffusion from the ISPH runs compared to the GDSPH runs, even with the highest resolution. There also seem to be additional noise in the small-scale density structures as we increase the number of neighbors for ISPH. Potentially related to the issue found in Fig. A.3

In the text

Fig. A.3 Particle density rendering of Kelvin-Helmholtz instability swirl at t = 2 for ∆ρ = 8. The figures shows the different behavior of not using diffusion versus using diffusion for the GDSPH and ISPH methods. A significant number of entropy clumps can be seen being generated in the low density medium for ISPH, as the removal of linear gradient together with sharp boundaries causes chaotic noise at the boundary. Lagrangian SPH formalisms that are derived directly from the density estimate, provide full spatial conservation of entropy. However, deviation from the Lagrangian formalism, while still using the traditional density estimate means that we will introduce entropy errors in our solution. GDSPH corrects for firstorder errors in entropy making it second-order in entropy. ISPH removes linear-gradient errors and will thus only be first-order accurate with entropy. The issue might be further aggravated by the use of average gradient kernels within ISPH. Adding thermal diffusion allows for local mixing between the cold and hot phase, and we can see that we get good behavior in both ISPH and GDSPH. Slight numerical surface tension effect can be seen for GDSPH, while none can be seen for ISPH.

In the text

In the text

In the text