Local Consistency of Smoothed Particle Hydrodynamics (SPH) in the Context of Measure Theory

The local consistency of the method of Smoothed Particle Hydrodynamics (SPH) is proved for a multidimensional continuous mechanical system in the context of measure theory. The Wasserstein distance of the corresponding measure-valued evolutions is used to show that full convergence is achieved in the joint limit N → ∞ and h → 0, where N is the total number of particles that discretize the computational domain and h is the smoothing length. Using an initial local discrete measure given by ${\mu }_0^N={\sum }_{b=1}^{N}m({\text{x}}_b,h){\delta }_{0,{\text{x}}_b(0)}$, where m_b = m(x_b, h) is the mass of particle with label b at position x_b(t) and δ_0,xb(t) is the x_b(t)-centered Dirac delta distribution, full consistency of the SPH method is demonstrated in the above joint limit if the additional limit ${N}$ → ∞ is also ensured, where ${N}$ is the number of neighbors per particle within the compact support of the interpolating kernel.

1. Introduction

In the late 1970s, the method of Smoothed Particle Hydrodynamics (SPH) was introduced independently by Gingold and Monaghan [1] and Lucy [2] as a tool for the simulation of astrophysical flows. SPH is a Lagrangian method for solving the equations of hydrodynamics, where the continuum is represented by a finite collection of point masses, or particles, which, in addition to mass, carry information of any other properties of the physical system. Since 1993, after being widely used in different applications in astrophysics, the method has found applications in a ever increasing number of areas in science and engineering. As an example, it has been used in such dissimilar problems as the friction dynamics of massive black holes [3] and the Taylor-Green vortex problem [4], just to mention a few. Even though the method is always gaining more practitioners, there still remain important drawbacks related to its poor convergence properties and loss of mathematical consistency. However, in the course of the years a notable amount of work has been devoted to solve the problem of SPH consistency [5–14]. In particular, the loss of consistency in standard SPH occurs because of zeroth-order truncation errors that arise when passing from the kernel to the particle approximation [12, 14]. Different sources of particle inconsistency have been previously identified by Liu et al. [15] to be associated with irregular particle distributions, kernel deficiencies at physical boundaries, and variable smoothing lengths. A new source related to the traditionally low number of neighbors within the kernel support was more recently discovered by Zhu et al. [12]. An alternative methodology that decreases the SPH discretization error to machine precision was introduced by da Silva et al. [16]. The method uses a post-processing technique based on repeated Richardson extrapolation (RRE) to improve the accuracy of the numerically obtained SPH solution. With an almost null cost in CPU time and memory usage, they obtained truncation error magnitudes of the same level of the machine round-off error for both a steady and unsteady heat diffusion model. Such level of error implies that consistency is actually restored to machine precision.

Measure theory has provided a powerful mathematical framework to successfully address the convergence of SPH by estimating the distance between a discrete measure and a smoothed measure of the continuum, when the total number of particles approaches infinity [5, 6, 13]. An early attempt to prove the convergence of SPH was reported by Di Lisio et al. [6], who used measures in combination with the Wasserstein distance. In particular, they employed the Wasserstein distance between a discrete measure and a regularized measure adopting a kernel function of extended domain and developed a proof of convergence for the regularized version of Euler's equations for a polytropic fluid model. Recently, Evers et al. [13] generalized the formalism of Di Lisio et al. and obtained a convergence proof for standard SPH [17]. They provided the order of convergence of SPH for a general class of force fields, including external and internal conservative forces, frictional dissipation, and non-local interactions. In this paper, we extend the convergence theorems derived by Evers et al. [13] to interpolation kernel functions of local finite (i.e., compact) support. We demonstrate that it is possible to limit the Wassertein distance between the initial discrete measure, associated with a non-uniform distribution of particles, and a continuous measure, associated with a regularization of a kernel with compact support.

The paper is organized as follows. The fundamentals of the SPH method are briefly described in Section 2. The convergence theorems of SPH in the context of measure theory are derived in Section 3. The Wasserstein distance between the discrete and the continuous measures is calculated for a non-uniform particle distribution and the convergence theorems proved by Evers et al. [13] are extended to interpolation kernel functions of compact support. The main conclusions are summarized in Section 4.

2. The SPH Method

2.1. Action for a Continuous Setting

In the standard SPH formalism [17], the derivation of the SPH equations of motion for a system of point masses starts from the discrete version of the Lagrangian L(t) for a continuous setting [13, 18], namely

$\begin{array}{ll}L(t):={\displaystyle {\int }_{{\Omega }_{s,t}}}\left[\frac{1}{2}\text{u}·\text{u}-e(\rho (\text{y},t),\text{y})\right]\rho (\text{y},t)d^n\text{y},& (1)\end{array}$

where u is the velocity field and ρ is the density, both defined as a function of time, t ∈ ℝ⁺ ≡ [0, ∞), and position, y ∈ ℝⁿ, where ℝⁿ denotes an n-dimensional Eulerian space, while e is the internal energy, defined as a function of density and position. The integration is performed over a volume domain ${\Omega }_{s,t}\subset {ℝ}^n$, which may deform during the time evolution of the system. As the integration is carried over space, the Lagrangian is only a function of time.

The trajectory of a material particle of the system, initially located at x = Φ(x, 0), can be defined by the coordinate transformation y = Φ(x, t) so that

$\begin{array}{ll}\text{u}(\Phi (\text{x},t),t)=\frac{\text{d}}{\text{d}t}\Phi (\text{x},t).& (2)\end{array}$

According to this motion mapping, the integral above can be transformed into

$\begin{array}{l}L\left[\Phi \right](t):\\ ={\displaystyle {\int }_{{\Omega }_{s,0}}}\left[\frac{1}{2}|\frac{\text{d}}{\text{d}t}\Phi (\text{x},t){|}^2-e(\rho (\Phi (\text{x},t),t),\Phi (\text{x},t))\right]\rho (\text{x},0)d^n\text{x},& (3)\end{array}$

where the coordinate transformation Φ(x, t) is such that Ω_s,t = Φ(Ω_s,0), where Ω_s,0 denotes the initial domain. In writing Equation (3), it has been assumed that ρ(x, t) and ρ(x, 0), as defined in the initial domain Ω(s, 0), are related by the same transformation mapping Φ(x, t). Invoking mass conservation, this relation reads as follows [19].

$\begin{array}{ll}\rho (\text{x},0)=\rho (\Phi (\text{x},t),t)|J\Phi (\text{x},t)|,& (4)\end{array}$

where |JΦ(x, t)| denotes the determinant of the Jacobian matrix of the coordinate transformation (2) and the integration is performed over the entire initial volume Ω_s,0 occupied by the continuous system. The action of the continuous system is given by the integral

$\begin{array}{ll}S\left[\Phi \right]={\displaystyle {\int }_0^T}L\left[\Phi \right](t)dt.& (5)\end{array}$

The equations of motion in SPH form can be derived by the principle of least action, that is, by minimizing the action (5) [17, 18]. The minimization of the action leads to the Euler-Lagrange equations, which for the particle discrete version of the Lagrangian, lead to the particle discretization of the equations of hydrodynamics. The smoothing and discretization procedure applied on a smooth function is called SPH function interpolation, or simply particle estimate of the function [14, 20]. The SPH discretization of the governing differential equations is performed in a two-step process. The first step is known as the kernel or smoothing approximation and the second step is the particle or SPH approximation. In what follows we shall provide a brief description of both approximations.

2.2. Kernel Approximation

The smoothing, or kernel approximation, consists of a weighted average of a field A(x) via an interpolation function or kernel W(|x − x′|, h)

$\begin{array}{ll}Ã(\text{x})={\displaystyle {\int }_{{\Omega }_n}}A({\text{x}}^{\prime })W(|\text{x}-{\text{x}}^{\prime }|,h)d^n{\text{x}}^{\prime },& (6)\end{array}$

where Ã(x) is the kernel estimate of the field A(x) and ${\Omega }_n\subseteq {ℝ}^n$ is the volume of the integration domain. The parameter h is called the smoothing length and determines the volume of influence of the smoothing function [21]. It is operationally defined as twice the standard deviation of the kernel [22], i.e.,

$\begin{array}{ll}h=2n^{1/2}{\left[{\displaystyle {\int }_{{ℝ}^n}}(\text{x}·\text{x})W(\left|\text{x}\right|,h)d^n\text{x}\right]}^{1/2},& (7)\end{array}$

where as before n denotes the dimension of the Euclidean space.

The kernel function W must satisfy the following conditions [14, 18, 20, 22, 23]:

(i) It must be positive definite, symmetric, and depend on the smoothing kernel and the distance between the field point and its neighbors.

(ii) In the limit h → 0, the kernel must tend to the Dirac delta distribution.

(iii) The kernel must be doubly continuous and differentiable. Furthermore, for all h it takes its maximum value at the location of the field point (x = 0) and must be a monotonically decreasing function everywhere.

(iv) It must satisfy the normalization condition

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_n}}W(|\text{x}-{\text{x}}^{\prime }|,h)d^n{\text{x}}^{\prime }=1.\end{array}$

(v) In modern applications most kernels have compact support, ${\Omega }_{n,\text{x}}\subset {ℝ}^n$ centered at x, so that only the attributes of neighboring particles contained within the sphere of influence of radius H and volume ${V}$_n,x are taken into account in the interpolation procedure. Hence W = 0 if |x − x′| ≥ H, where H = kh and k is a constant factor that determines the size of the kernel support.

2.3. Particle Approximation

In the particle approximation the domain is divided into N Lagrangian subdomains, or cells, in much the same way as for traditional grid-based methods. A particle is placed at the center of each cell and the boundaries of the cells are such that the particle masses remain always constant. Using the mean value theorem, the particle estimate of the smoothed field Ã(x) defined by Equation (6) at the position x_a of particle a is given by the approximation

$\begin{array}{ll}{A}_a={\displaystyle \sum _{b=1}^{{N}({\text{x}}_a,h)}}{A}_{b}W(|{\text{x}}_a-{\text{x}}_b|,h)\Delta V_b,& (8)\end{array}$

where x_b (for b = 1⋯${N}$(x_a, h)) denotes the positions of neighboring particles within the compact support ${\Omega }_{n,{\text{x}}_a}\subset {ℝ}^n$ of particle a, ΔV_b is the volume associated to particle b, and ${N}$(x_a, h) is the average number of particles enclosed by the compact support of the kernel. It is common practice in SPH to replace the volume ΔV_b by the ratio m_b/ρ_b, where m_b is the mass of particle b and ρ_b is its density. For a non-uniform distribution of particles ${N}$(x_a, h) can be expressed as [14],

$\begin{array}{ll}{N}({\text{x}}_a,h)={{V}}_{n,a}\frac{\rho ({\text{x}}_a)}{m({\text{x}}_a,h)}+O\left(h^{n+2}\right),& (9)\end{array}$

when h ≪ 1. Kernels with a compact support are characterized by two parameters: h and H = kh, where H is the effective radius of the spherical support [14, 22]. Therefore, the dependence of ${N}$(x_a, h) on kh is evident. The dominant term in Equation (9) is just the quotient between the volume of the compact support, ${V}$_n,a, and the volume ΔV_a of particle a. This ratio defines the total number of neighbors of particle a up to order n + 2 in the smoothing length. For sufficiently small values of h this provides a very good approximation to determine the average number of neighbors per particle.

If the positions x_b (with b = 1, 2, …, ${N}$(x_a, h)) of all neighbors of particle a are known then we can obtain maximum information on their number and distribution. For example, the distance between pairs of particles, Δ(x_a, x_b), which is a scalar correlation measure of how the particles are distributed within the compact support, is a set of ${N}$(x_a, h)[${N}$(x_a, h)−1]/2 elements over which we can construct as measures the minimum distance, Δ_min, the maximum distance, Δ_max, and the average distance, Δ_mean, between all countable pairs of neighbors within the compact support of the kernel. The mean distance Δ_mean between pairs of particles can be defined as

$\begin{array}{ll}{\Delta }_{mean}={\left[\frac{{{V}}_{n,a}}{{N}({\text{x}}_a,h)}\right]}^{1/n}.& (10)\end{array}$

The position x_b of particle identified by label b can be interpreted as a mapping between a Euclidean space ℝⁿ of configurations, where the particles take their possible positions, and a vector label space b ∈ ℝⁿ. The determinant of the Jacobian matrix J_xb(b) of the mapping x_b in the context of the SPH interpolation is given by [14].

$\begin{array}{ll}\left|J_{{\text{x}}_b}(\text{b})\right|=\left|\frac{\partial (x_{1b},x_{2b},\cdots ,x_{nb})}{\partial (b_1,b_2,\cdots ,b_n)}\right|=\frac{m_b}{{\rho }_b}\ne 0,& (11)\end{array}$

where b = b(x_b) is a bijective map that admits the inverse x_b = x_b(b) and the kernel of the transformation is of zero dimension. Therefore, two particles x_a and x_b, with labels a ≠ b, are two disjoint particles without interstices between them (see Figure 1). If kernel functions of compact support are used to perform the SPH interpolation, Sigalotti et al. [14] demonstrated that partition of unity can be achieved only up to the order O(hⁿ⁺²), i.e.,

$\begin{array}{ll}\frac{1}{{{V}}_{n,a}}\left({\displaystyle \sum _{b=1}^{{N}\left({\text{x}}_a,h\right)}}\Delta V_b\right)=1+O\left(h^{n+2}\right).& (12)\end{array}$

That is, the ratio of the sum of the volumes of all neighbors over the volume of the kernel support is exactly unity only when h → 0 and ${N}$(x_a, h) → ∞, in which case N → ∞ also. This can be easily observed from Figure 1, where neighbors located close to the border of the compact support may contribute by excess or defect to the sum of the particle volumes ΔV_b.

FIGURE 1

Figure 1. Schematic drawing illustrating the SPH interpolation in two-space dimensions (n = 2). The circle of radius H = kh defines the compact support of the kernel function W(|x_a − x_b|, h) centered at the position x_a of the observation particle a. Although the particles are depicted as point-like masses, their geometrical shape corresponds to an irregular rectangular cell of area m_b/ρ_b. Due to the Lagrangian nature of SPH, the rectangles deform during the evolution. However, they do not overlap and do not have empty interstices between them, which is consistent with the requirement that m_b/ρ_b ≠ 0.

2.4. Scaling Laws and Resolution of the SPH Interpolation

In most SPH simulations it has been a common practice to fix the resolution by simply choosing the value of h as a factor of the interparticle distance. However, Zhu et al. [12] derived scaling relations complying with the resolution limit h → 0 when N → ∞. If a kernel of compact support is employed, then full consistency for the SPH interpolation is achieved when the limit ${N}$(x_a, h) → ∞ is also satisfied. Only when this joint limit is met the particle approximation, defined by Equation (8), converges locally to the kernel approximation, defined by Equation (6).

According to Equation (9) and for a non-uniform distribution of particles, the limit ${N}$(x_a, h) → ∞ is achieved for h → 0 if the mass, m(x_a, h), of particles scales with h as h^β with β > n [14]. Consequently, a local scaling law is obtained for the number of neighbors ${N}$(x_a, h) within the compact support Ω_n,a of the kernel as

$\begin{array}{ll}{N}({\text{x}}_a,h)~h^{n-\beta },\text{ for }\beta >n.& (13)\end{array}$

This is the same scaling law obtained by Zhu et al. [12] for the case when n = 3. A global scaling law for the total number of particles N is obtained by invoking a relation of similarity between the average particle density in the compact support Ω_n,a of the kernel and the particle density in the full computational domain of total volume V:

$\begin{array}{ll}\frac{{N}({\text{x}}_a,h)}{{{V}}_{n,a}}=\frac{N}{V}.& (14)\end{array}$

Since the volume V of the computational domain does not depend on the parameters of the SPH interpolation, the global scaling law is obtained as follows

$\begin{array}{ll}N~h^{-\beta },\text{ for }\beta >n,& (15)\end{array}$

which complies with the joint limit for full particle consistency.

3. SPH Convergence Theorem

A measure-value formulation was used by Evers et al. [13] to prove the convergence of SPH as a special case of a class of approximating measures that is much broader than just a sum of weighting kernel functions. In their case the proof was obtained in the limit when N → ∞ for kernel functions of extended support. Here we generalize the measure-value formulation implemented by Evers et al. [13] to prove the theorem of SPH convergence for kernel functions of locally finite support ${\Omega }_{n,\text{y}}\subset {ℝ}^n$, where y is the coordinate in a Eulerian description. The analysis applies equally well to regularly and irregularly distributed particles. For a fixed instant of time t ∈ [0, T], where T defines the evolution time, the compact support of the kernel will be denoted by Ω_y,t. Motion mapping, y = Φ(x, t), is invoked (with the initial condition x = Φ(x, 0)), in order to transform from Eulerian (y) to Lagrangian (x) coordinates (see Equation 2).

3.1. SPH Interpolation in the Context of Measure Theory

In the context of measure theory, the mass differential ρ(y, t)dⁿy is associated with the measure ${\mu }_t(d^n\text{y})$ for a given time t, while at time t = 0 the mass differential ρ(x, 0)dⁿx is associated to the measure ${\mu }_0(d^n\text{x})$. Therefore, the kernel approximation for the mass density of the fluid is

$\begin{array}{ll}{\tilde{\rho }}_h\left(\text{y},t\right)={\displaystyle {\int }_{{\Omega }_{\text{y},t}}}W\left(|\text{y}-{\text{y}}^{\prime }|,h\right){\mu }_t\left(d^n\text{y}\right).& (16)\end{array}$

In this regularization process the kernel function W(|y − y′|, h) has compact support Ω_y,t and meets the conditions listed in Section 2.2. Note that the regularized density depends on the smoothing length h. If we invoke the motion mapping y = Φ(x, t) (or motion push-forward) in the context of measure theory, then the kernel approximation of the mass density ${\tilde{\rho }}_h(\text{y},t)$ becomes

$\begin{array}{ll}{\tilde{\rho }}_h\left(\Phi \left(\text{x},t\right),t\right)={\displaystyle {\int }_{{\Omega }_{\text{x},0}}}\left(|\Phi \left(\text{x},t\right)-\Phi \left({\text{x}}^{\prime },t\right)|,h\right){\mu }_0\left(d^n{\text{x}}^{\prime }\right).& (17)\end{array}$

Equation (16) for ${\tilde{\rho }}_h(\text{y},t)$ is a functional of μ_t(•), while Equation (17) for ${\tilde{\rho }}_h(\Phi (\text{x},t),t)$ is a functional of Φ(x, t).

The particle approximation is obtained by substituting the measure μ₀(•) by the discrete measure ${\mu }_0^N(•)$ so that

$\begin{array}{ll}{\mu }_0^N={\displaystyle \sum _{b=1}^N}{m}_b{\delta }_{{\text{x}}_b,0},& (18)\end{array}$

where the mass m_b of particle b scales with h as h^β, with β > n [14]. The discrete measure ${\mu }_t^N(•)$ is constructed using the Dirac delta measures δ_0,xb(t) centered at x_b(t), which evolve according to the push-forward motion, x_b(t) = Φ(x_b, t), with x_b(0) = x_b. Substitution of the initial discrete measure ${\mu }_0^N(•)$ given by Equation (18) into the Lagrangian L(t) of the continuous system given by Equation (3) and in the regularized density ${\tilde{\rho }}_h(\Phi (\text{x},t),t)$ defined by Equation (17), we obtain the discrete version of the Lagrangian for a system of point masses as it applies to the traditional SPH method [17, 22], i.e.,

$\begin{array}{ll}L(t)={\displaystyle \sum _{b=1}^N}{m}_b\left[\frac{1}{2}|{\dot{\text{x}}}_b(t){|}^2-e({\tilde{\rho }}_h({\text{x}}_b(t),t),{\text{x}}_b(t))\right],& (19)\end{array}$

with

$\begin{array}{l}{\tilde{\rho }}_h({\text{x}}_a(t))={\displaystyle \sum _{b=1}^{{N}({\text{x}}_a(t),h)}}{m}_{b}W(|{\text{x}}_a(t)-{\text{x}}_b(t)|,h),\end{array}$

where N is the total number of particles that discretize the computational volume and ${N}$(x_a(t), h) is the number of neighbors within the compact support of the kernel W(|x_a(t) − x_b(t)|, h) centered at the position x_a(t) of the observation particle a.

3.2. The Limit ${\mu }_0^N\to {\mu }_0$ When N → ∞

In order to compute the limit of the discrete measure ${\mu }_0^N$ when N → ∞ it is necessary to define eigenobjects in the measure theory, such as the push-forward, the joint representation, and the Wasserstein distance [13].

Definition 3.1 (Push-forward). For a map Φ:A₁ → A₂ between two measurable spaces (A₁, σ₁) and (A₂, σ₂) the push-forward for the measure $\mu :{\sigma }_1\to 0\cup {ℝ}^{+}$ is the measure that is constructed by means of the map Φ:

$\begin{array}{l}\left[\Phi \#\mu \right](B)=\mu \left[{\Phi }^{-1}(B)\right]\forall B\in {\sigma }_2,\end{array}$

or in its integral form

$\begin{array}{l}{\displaystyle {\int }_B}f(\text{y})\left[\Phi \#\mu \right](d^n\text{y})={\displaystyle {\int }_{{\Phi }^{-1}(B)}}f(\Phi (\text{x}))\mu (d^n\text{x}),\end{array}$

where the map Φ defines the change of coordinates y = Φ(x, t).

The push-forward definition contains the change of variables in the integration methods, as observed in the change from the Eulerian to the Lagrangian representation through the motion mapping y = Φ(x, t) (see Equations 1–3, 16, and 17). In particular, the particle approximation can be expressed as

$\begin{array}{ll}{\Phi }_{\text{SphPA}}:V\subset {ℝ}^n\to \left\{\Delta V_b,\text{with }b=1,2,\dots ,N\right\},& (\mathrm{20a})\end{array}$

with

$\begin{array}{ll}{\Phi }_{\text{SphPA}}(\text{x})={\text{x}}_b\text{ if }\text{x}\in \Delta V_b\text{ for some }b\in \left\{1,2,\dots ,N\right\},& (\mathrm{20b})\end{array}$

where ΔV_b is the volume of particle b and x_b is its position. Then

$\begin{array}{l}{\displaystyle {\int }_V}f(\text{x})\left[{\Phi }_{\text{SphPA}}\#{\mu }_0\right](d^n\text{x})={\displaystyle {\int }_V}f({\Phi }_{\text{SphPA}}(\text{x})){\mu }_0(d^n\text{x}),\end{array}$

with $\dots {\mu }_0^n(d^n\text{x})=\dots \left[{\Phi }_{\text{SphPA}}\#{\mu }_0\right](d^n\text{x})=\sum _{b=1}^N{m}_b{\delta }_{{\text{x}}_b,0}\dots d^n\text{x}$.

Note that for all x ∈ V there always exists b that belongs to the set of labels {1, 2, …, N} such that

$\begin{array}{ll}|\text{x}-{\Phi }_{\text{SphPA}}(\text{x})|=|\text{x}-{\text{x}}_b|<\alpha {\Delta }_{mean}.& (21)\end{array}$

For bounded V, there exists α ∈ ℝ⁺ such that Δ_sup = αΔ_mean, where Δ_sup is the maximum distance between all pairs of particles Δ(x_a, x_b) for a, b = 1, 2, …, N. In the present analysis we consider the parameter α of the distribution of particles and not of the total number of particles N. This is consistent with the similarity relationship (14).

The partition of unity as introduced by Sigalotti et al. [14] can be expressed in the context of measure theory. Let ${V}$_n,a be the spherical volume of the compact support Ω_n,a of the kernel function $W(|{\text{x}}_a-{\text{x}}^{\prime }|,h)$, centered at position x_a and with smoothing length h, then the partition of unity can be written as

$\begin{array}{ll}\frac{{\displaystyle {\int }_{{{V}}_{n,a}}}{d}^n\text{x}}{{{V}}_{n,a}}=\frac{{\displaystyle {\int }_{{V}n,a}}\left[\frac{1}{{\rho }_0(\text{x})}\right]{\mu }_0(d^n\text{x})}{{{V}}_{n,a}}=1,& (22)\end{array}$

where the domain of integration is taken over the volume of the kernel support. When passing from the smoothed to the particle approximation, the partition of unity will then take the form [14]:

$\begin{array}{l}\frac{{\displaystyle {\int }_{{V}n,a}}\left[\frac{1}{{\rho }_0({\Phi }_{\text{SphPA}}(\text{x}))}\right]{\mu }_0(d^n\text{x})}{{{V}}_{n,a}}\to \frac{{\displaystyle {\sum }_{b=1}^{{N}({\text{x}}_a,h)}}\Delta V_b}{{{V}}_{n,a}}=1\end{array}$

$\begin{array}{ll}+O\left(h^{n+2}\right),& (23)\end{array}$

where ${N}$(x_a, h) is the number of neighbors that belong to Ω_n,a and ΔV_b = m_b/ρ(x_b) is the volume of particle b located at the point x_b, with x_b ∈ Ω_n,a. In the particle approximation the partition of unity is reproduced only approximately to order O(hⁿ⁺²). In Figure 1, where the kernel support is represented schematically in the plane (n = 2), the partition of unity deviates from one due to those particle volumes that are intersected by the circular compact support.

Definition 3.2 (Joint representation). Let two measures be ${\mu }_1:\sigma \to 0\cup {ℝ}^{+}$ and ${\mu }_2:\sigma \to 0\cup {ℝ}^{+}$ of a space of measure (A, σ). The joint representation of μ₁ and μ₂ is the measure π:σ × σ → 0∪ℝ⁺ such that

$\begin{array}{l}{\displaystyle {\int }_{\sigma \times \sigma }}f(x_i)\pi (d{x}_1,d{x}_2)={\displaystyle {\int }_{\sigma }}f(\text{x}){\mu }_i(dx),i=1,2.\end{array}$

Definition 3.3 (Wasserstein distance). Let Π(μ₁, μ₂) be the set of all joint representations of μ₁ and μ₂. The Wasserstein distance between two measures μ₁ and μ₂ of a measurable space (A, σ) is

$\begin{array}{l}\begin{array}{c}{W}({\mu }_1,{\mu }_2):={\displaystyle \underset{\pi \in \Pi ({\mu }_1,{\mu }_2)}{\text{inf}}}{\displaystyle {\int }_{\sigma \times \sigma }}|x-y|\pi (dx,dy).\end{array}\end{array}$

Using Equation (21), the Wasserstein distance between the discrete measure ${\mu }_0^N$ and the continuous measure μ₀ is

$\begin{array}{l}{W}({\mu }_0^N,{\mu }_0)\le {\displaystyle {\int }_V}|\text{x}-{\Phi }_{\text{SphPA}}(\text{x})|{\mu }_0(d^n\text{x})<\end{array}$

$\begin{array}{ll}\alpha {\Delta }_{mean}{\displaystyle {\int }_V}{\mu }_0(d^n\text{x})<\alpha V{\Delta }_{mean}.& (24)\end{array}$

In the limit N → ∞, the result of Evers et al. [13] is obtained from the above inequality provided that the similarity relation (14) holds. Hence

$\begin{array}{ll}{W}({\mu }_0^N,{\mu }_0)\le {\displaystyle {\int }_V}|\text{x}-{\Phi }_{\text{SphPA}}(\text{x})|{\mu }_0(d^n\text{x})<\alpha {\left[\frac{V^{n-1}}{N}\right]}^{1/n},& (25)\end{array}$

and therefore taking the limit when N → ∞ gives

$\begin{array}{ll}{\displaystyle {lim}_{N\to \infty }}{W}({\mu }_0^N,{\mu }_0)\to 0.& (26)\end{array}$

Furthermore, if we add the necessary requirement to satisfy the global scaling (15), we immediately find that

$\begin{array}{ll}{W}\left({\mu }_0^N,{\mu }_0\right)={W}\left({\mu }_0^h,{\mu }_0\right)<h^{\beta /n}\text{with}\beta >n,& (27)\end{array}$

so that in the limit when h → 0

$\begin{array}{ll}{\displaystyle {lim}_{h\to 0}}{W}\left({\mu }_0^h,{\mu }_0\right)\to 0.& (28)\end{array}$

3.3. Dynamics of the Regularized Continuous System

According to Evers et al. [13] the Lagrangian L(t) (see Equation 3) is regularized as follows. First, ρ(x, 0)dⁿx is replaced by μ(dⁿx) and then the density is regularized according to Equation (17) for ${\tilde{\rho }}_h(\Phi (\text{x},t),t)$, where now the kernel function W(|Φ(x, t) − Φ(x′, t)|, h) has compact support initially centered at x.

The dynamics of the regularized system is obtained by generalizing the principle of least action [13]. The dynamics obtained this way describe the evolution of two fields: the trajectory of the material point Φ(x, t) and the measure μ_t. The dynamic state is well-defined if the initial position of each material point, namely Φ(x, 0) = x, its initial velocity field ${\stackrel{•}{\Phi }}_0=v_0(\text{x})$, and the initial measure μ₀, which is associated with the geometry of the space where the continuous system evolves, are all known.

Then the dynamical system to be described is:

Given T > 0, we have a V ⊂ ℝⁿ and a σ-algebra, denoted by Σ, of the set V with the initial measure μ₀ over the Measurable Space (V, Σ) that satisfies μ₀ : Σ → [0, +∞) and an initial velocity field $v_0\in C^1({ℝ}^n;{ℝ}^n)$. The evolution is then governed by the following set of equations

$\begin{array}{ll}\{\begin{array}{l}\ddot{\Phi }(x,t)=-\frac{\partial e}{\partial \tilde{\rho }}\left[{\tilde{\rho }}_h\left(\Phi (x,t)\right)\right]\nabla {\tilde{\rho }}_h\left(\Phi (x,t)\right)-\left(\nabla W_h\ast \left[\left(\frac{\partial e}{\partial \tilde{\rho }}°{\tilde{\rho }}_h\right){\mu }_t\right]\right)\left(\Phi (x,t)\right)\hfill \\ -\nabla U\left(\Phi (x,t)\right)-\eta \left(\Phi (x,t)\stackrel{•}{\Phi }(x,t)\right)+\left(K\ast {\mu }_h\right)\left(\Phi (x,t)\right),\hfill \\ {\tilde{\rho }}_h:=W_h\ast {\mu }_t,\hfill \\ {\mu }_t=\Phi (x,t)\#{\mu }_0,\hfill \\ \Phi (x,0)=x;\dot{\Phi }(x,0)=v_0(x).\hfill \end{array}& (29)\end{array}$

For all x ∈ V and all t ∈ (0, T] the notation W_h * μ_t emphasizes that the evolution in the regularization is caused by the push-forward Φ(x, t) # μ₀, where U ∈ C²(ℝⁿ; ℝ) is the potential due to external fields. The internal energy e(ρ) ∈ C²(ℝ⁺; ℝ) is associated only with the thermodynamic properties of the continuum, with e(ρ, y) = U(y)+e(ρ). The field η ∈ C¹(ℝⁿ ℝ⁺) describes the viscosity and the field K ∈ C¹(ℝⁿ ℝⁿ) represents a non-local internal interaction in the continuous system.

3.4. The Theorem of Existence-Uniqueness and Continuity With Respect to the Initial Conditions

The main theorem of Evers et al. [13] ensures that the regularized dynamical system given by Equation (29), with the initial conditions of Φ(x, 0) = x and $\dot{\Phi }(\text{x},0)=v_0(\text{x})$, has a unique solution (μ_t, Φ(x, t)) for all x ∈ V and t ∈ (0, T], if the initial measure μ₀ is known. Also the solutions of Equation (29) are continuous with respect to the initial measure. The following theorems hold:

Theorem 3.1 (Existence-uniqueness). For all t ∈ (0, T], let h > 0 and the initial conditions Φ(x, 0) = x and $\dot{\Phi }(\text{x},0)=v_0(\text{x})$ known for the regularized continuous system of Equation (29). Then for each initial measure μ₀ over the measurable space (V, Σ) there is a unique pair solution (μ_t, Φ(x, t)) to the dynamical system defined by Equation (29).

Theorem 3.2 (Continuity with respect to μ₀). For all t ∈ (0, T], let be h > 0 and the initial conditions Φ(x, 0) = x and $\dot{\Phi }(0,\text{x})=v_0(\text{x})$ known for the regularized continuous system of Equation (29). Given two initial measures ${\mu }_0^1$ and ${\mu }_0^2$ over the measurable space (V, Σ), the supremum of the Wasserstein distance between the two solutions $({\mu }_t^1,{\Phi }^1(\text{x},t))$ and $({\mu }_t^2,{\Phi }^2(\text{x},t))$ of the regularized dynamical system given by Equation (29) is bounded by:

$\begin{array}{l}{\displaystyle \underset{t\in \left[0,T\right]}{\text{sup}}}{W}\left({\mu }_t^1,{\mu }_t^2\right)\le \left(1+T\left(\Vert \nabla v_0{\Vert }_{\infty }+\Vert \eta {\Vert }_{\infty }\right)\right)\\ exp\left(\Vert \eta {\Vert }_{\infty }T+\frac{1}{2}\left(M_4+M_5\right)T^2\right){W}\left({\mu }_0^1,{\mu }_0^2\right),\end{array}$

where M₄ and M₅ are both positive and independent of the initial measurements ${\mu }_0^1$ and ${\mu }_0^2$.

Corollary 3.1 (Convergence of the particle approximation to the smoothing approximation). For all t ∈ (0, T], let be h > 0 and the initial conditions be Φ(x, 0) = x and $\dot{\Phi }(0,\text{x})=v_0(\text{x})$ known for the regularized continuous system of Equation (29). Consider a sequence of initial discrete measures ${\left\{{\mu }_0^N\right\}}_{N\in \mathbb{N}}$ over the measurable space (V, Σ) that converge to measure μ₀ over (V, Σ), then

$\begin{array}{l}{W}({\mu }_0^N,{\mu }_0){\displaystyle \underset{}{\overset{N\to \infty }{\to }}}0.\end{array}$

For each sequence of initial discrete measures, $({\mu }_0^N,\Phi (\text{x},0))$, there is a unique sequence $({\mu }_t^N,{\Phi }^N(\text{x},t))$ solution of the dynamical system given by Equation (29). Let (μ_t, Φ(x, t)) be the solution of the dynamical system of Equation (29) with initial conditions (μ₀, Φ(x, 0)). Then ${\mu }_t^N$ converges to μ_t as

$\begin{array}{ll}{\displaystyle \underset{t\in \left[0,T\right]}{\text{sup}}}{W}({\mu }_t^N,{\mu }_t){\displaystyle \underset{}{\overset{N\to \infty }{\to }}}0.& (30)\end{array}$

Remark (On the convergence of the particle approximation to the kernel approximation when h → 0). Corollary 3.1 ensures that in the limit h → 0 the SPH particle approximation, Φ_SphPA (see Equation 20), converges to the kernel approximation. This is true if the local scaling law (13) and the similarity relation (14) are satisfied [or equivalently, the similarity relation (14) and the global scaling law (15) are both satisfied].

4. Analytical Proof of Convergence

Full consistency of the particle approximation is guaranteed only when the joint limit N → ∞, h → 0 and ${N}$ → ∞ is achieved [12]. Previous studies on the SPH consistency were mostly based on how well the discrete kernel normalization condition is reproduced [12, 14]. For low discrepancy sequences of particles, the standard deviation of the density field approximately varies as ${N}$⁻¹ when ${N}$ is larger than about 200 [12]. According to Corollary 3.1, the Wasserstein distance can be used to measure how the particle approximation converges to the continuous kernel approximation. Let us first consider the case of an n-dimensional system of volume V, for which the kernel approximation is implemented using a non-local kernel function (i.e., without compact support). The Wasserstein distance between a discrete measure ${\mu }_t^N$ and a continuous measure μ_t can be estimated from its upper bound given by Equation (25) as

$\begin{array}{ll}{W}({\mu }_t^N,{\mu }_t)\approx \alpha {\left[\frac{V^{n-1}}{N}\right]}^{1/n}.& (31)\end{array}$

If, on the other hand, the kernel approximation is implemented using an interpolation function of compact support, from Equation (27) it follows that

$\begin{array}{ll}{W}\left({\mu }_t^h,{\mu }_t\right)\approx h^{\beta /n}\text{with}\beta >n.& (32)\end{array}$

From the triangular inequality, the convergence properties of ${W}({\mu }_t^N,{\mu }_t)$, with t ∈ [0, T], can be studied by assessing that ${\mu }_t^N$ is a Cauchy sequence, that is

$\begin{array}{l}{W}\left({\mu }_t^{N_1},{\mu }_t^{N_2}\right)\le {W}\left({\mu }_t^{N_1},{\mu }_t\right)+{W}\left({\mu }_t^{N_2},{\mu }_t\right),\end{array}$

where N₂ > N₁ → ∞. Numerical calculations to assess that ${\mu }_t^N$ is a Cauchy sequence under the Wasserstein distance must require using finite linear programming techniques. However, the numerical calculation of the Wasserstein distance is a rather hard task because the computational complexity scales poorly as N. According to Panaretos and Zemel [24], the best algorithms to compute the optimal solution have the prohibitive complexity N³logN. Therefore, a proper numerical three-dimensional test showing the SPH convergence in terms of the Wasserstein distance for high values of N and ${N}$ and sufficiently low values of h will be the subject of a further research paper in this line.

When the SPH interpolation is performed using a kernel function of compact support, the Wasserstein distance between discrete measures, namely ${\mu }_t^{2^k}$ and ${\mu }_t^{2^{k+1}}$, for k ≫ 1 is given asymptotically by the relation

$\begin{array}{ll}{W}_{k,k+1}:=\frac{\left(\alpha V^{\frac{n-1}{n}}\right)\left(2^{\frac{1}{n}}+1\right)}{2^{\frac{k+1}{n}}}\ge {W}\left({\mu }_t^{2^k},{\mu }_t^{2^{k+1}}\right),& (33)\end{array}$

while the convergence rate for an n-dimensional system is given by

$\begin{array}{ll}{C}_{k+1}^n:={log}_2\left(\frac{W_{k+1,k+2}}{W_{k,k+1}}\right)=-\left(1/n\right),& (34)\end{array}$

where n is the dimension of the Euclidean space ℝⁿ. In particular, Evers et al. [13] have obtained numerically the convergence rate for n = 1 and 2, but using SPH parameters, i.e., low values of N and h, for which consistency cannot be guaranteed.

For a kernel function of a compact support, the Wasserstein distance between the two discrete measures ${\mu }_t^h$ and ${\mu }_t^{h+\Delta h}$ approaches asymptotically the form

$\begin{array}{ll}{W}_{h,h+\Delta h}:=h^{\beta /n}+{\left(h+\Delta h\right)}^{\beta /n}\ge {W}\left({\mu }_t^h,{\mu }_t^{h+\Delta h}\right),& (35)\end{array}$

for β > n when Δh ≪ h → 0. Hence, the rate of convergence obeys the relation

$\begin{array}{ll}\frac{d}{dh}{W}_{h,h+\Delta h}:=(\beta /n)h^{(\beta /n)-1}\left\{1+{(1+\Delta h/h)}^{(\beta /n)-1}\right\},& (36)\end{array}$

which embodies a rather slow convergence even when Δh/h < < 1 and h → 0.

5. Conclusions

In this paper the formulation introduced by Evers et al. [13] to prove the convergence of the method of Smoothed Particle Hydrodynamics (SPH) in the Wasserstein distance of the corresponding measure-valued evolutions is extended to kernel functions of locally finite support and irregular particle distributions. The present analysis has demonstrated the convergence and consistency of SPH in the joint limit h → 0 and N → ∞. In this limit consistency of the particle approximation solution of the equations of motion of the regularized continuum system is obtained when the number of neighbors within the compact support of the kernel, ${N}$, tends simultaneously to infinity, albeit at a rate slower than N so that ${N}$/N → 0. It is important to recall that, in contrast to previous proofs of SPH convergence in the framework of measure theory, here the continuous system considered is regularized by means of a kernel function of compact support. The particle approximation of the exact differential equations has been done in compliance with the similarity relation ${N}$/${V}$ = N/V, where ${V}$ is the volume enclosed by the kernel support and N is the total number of particles that discretize the computational domain, and the global scaling law N ~ h^−β for β > n, where h is the smoothing length and n is the dimension of the Euclidean space ℝⁿ. Convergence is also achieved if the above global scaling law is replaced by the local scaling ${N}$ ~ h^{n − β} for β > n, as was first derived by Zhu et al. [12] and later on by Sigalotti et al. [14] following different methodologies. Since most modern SPH simulations employ interpolation functions of compact support, the present analysis demonstrate that in the most favorable case when both N → ∞ and h → 0 simultaneously, then full consistency is achieved only when also ${N}$ → ∞.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

OR, GA, JK, LS, and CV have contributed intellectually to the form and content of this manuscript. The writing was done by OR and LS. All authors agree with the version submitted for publication, contributed to the article, and approved the submitted version.

Funding

We acknowledge funding from the European Union's Horizon 2020 Programme under the ENERXICO Project, grant agreement No. 828947 and from the Mexican CONACYT-SENER Hidrocarburos Programme under grant agreement B-S-69926.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We acknowledge financial support from the Departamento de Ciencias Básicas of the Universidad Autónoma Metropolitana - Azcapotzalco through internal funds.

References

1. Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Not R Astron Soc. (1977) 181:375–89. doi: 10.1093/mnras/181.3.375

CrossRef Full Text | Google Scholar

2. Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J. (1977) 82:1013–24. doi: 10.1086/112164

CrossRef Full Text | Google Scholar

3. Sayeb M, Blecha L, Kelley LZ, Gerosa D, Kesden M, Thomas J. Massive black hole binary inspiral and spin evolution in a cosmological framework. Monthly Not R Astron Soc. (2020) 501:2531–46. doi: 10.1093/mnras/staa3826

CrossRef Full Text | Google Scholar

4. Lind SJ, Stansby PK. High-order Eulerian incompressible smoothed particle hydrodynamics with transition to Lagrangian free-surface motion. J Comput Phys. (2016) 326:290–311. doi: 10.1016/j.jcp.2016.08.047

CrossRef Full Text | Google Scholar

5. Di Lisio R, Grenier E, Pulvirenti M. On the regularization of the pressure field in compressible Euler equations. Annali della Scuola Normale Superiore di Pisa, Classe di Scienze 4e. (1997) 24:227–38.

6. Di Lisio R, Grenier E, Pulvirenti M. The convergence of the SPH method. Comput Math Appl. (1998) 35:95–102. doi: 10.1016/S0898-1221(97)00260-5

CrossRef Full Text | Google Scholar

7. Ben Moussa B, Vila JP. Convergence of SPH method for scalar nonlinear conservation laws. SIAM J Num Anal. (2000) 37:863–87. doi: 10.1137/S0036142996307119

CrossRef Full Text | Google Scholar

8. Quinlan NJ, Basa M, Lastiwka M. Truncation error in mesh-free particle methods. Int J Num Methods Eng. (2006) 66:2064–85. doi: 10.1002/nme.1617

CrossRef Full Text | Google Scholar

9. Vaughan GL, Healy TR, Bryan KR, Sneyd AD, Gorman RM. Completeness, conservation and error in SPH for fluids. Int J Num Methods Fluids. (2008) 56:37–62. doi: 10.1002/fld.1530

CrossRef Full Text | Google Scholar

10. Read JI, Hayfield T, Agertz O. Resolving mixing in smoothed particle hydrodynamics. Monthly Not R Astron Soc. (2010) 405:1513–30. doi: 10.1111/j.1365-2966.2010.16577.x

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Fatehi R, Manzari MT. Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives. Comput Math Appl. (2011) 61:482–98. doi: 10.1016/j.camwa.2010.11.028

CrossRef Full Text | Google Scholar

12. Zhu Q, Hernquist L, Li Y. Numerical Convergence in smoothed particle hydrodynamics. Astrophys J. (2015) 800:6. doi: 10.1088/0004-637X/800/1/6

CrossRef Full Text | Google Scholar

13. Evers JHM, Zisis IA, van der Linden BJ, Duong MH. From continuum mechanics to SPH particle systems and back: Systematic derivation and convergence. J Appl Math Mech. (2018) 98:106–133. doi: 10.1002/zamm.201600077

CrossRef Full Text | Google Scholar

14. Sigalotti LDG, Rendon O, Klapp J, Vargas CA, Cruz F. A new insight into the consistency of the SPH interpolation formula. Appl Math Comput. (2019) 356:50–73. doi: 10.1016/j.amc.2019.03.018

CrossRef Full Text | Google Scholar

15. Liu MB, Liu GR. Restoring particle consistency in smoothed particle hydrodynamics. Appl Num Math. (2006) 56:19–36. doi: 10.1016/j.apnum.2005.02.012

CrossRef Full Text | Google Scholar

16. da Silva LP, Marchi CH, Maneguette M, Foltran AC. Robust RRE technique for increasing the order of accuracy of SPH numerical solutions. Math Comput Simul. (2022) 199:231–52. doi: 10.1016/j.matcom.2022.03.016

CrossRef Full Text | Google Scholar

17. Monaghan JJ. Smoothed particle hydrodynamics. Rep Prog Phys. (2005) 68:1703–59. doi: 10.1088/0034-4885/68/8/R01

CrossRef Full Text | Google Scholar

18. Price DJ. Smoothed particle hydrodynamics and magnetohydrodynamics. J Comput Phys. (2012) 231:759–94. doi: 10.1016/j.jcp.2010.12.011

CrossRef Full Text | Google Scholar

19. Landau LD, Lifshitz EM. Fluid Mechanics. Volume 6 of Course of Theoretical Physics. Oxford: Pergamon Press (1987).

20. Lind SJ, Rogers BD, Stansby PK. Review of smoothed particle hydrodynamics: towards converged Lagrangian flow modelling. Proc R Soc A Math Phys Eng Sci. (2020) 476:20190801. doi: 10.1098/rspa.2019.0801

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Liu GR, Liu MB. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. Singapore: World Scientific (2003) doi: 10.1142/9789812564405

CrossRef Full Text | Google Scholar

22. Dehnen W, Aly H. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Monthly Not R Astron Soc. (2012) 425:1068–82. doi: 10.1111/j.1365-2966.2012.21439.x

CrossRef Full Text | Google Scholar

23. Fulk DA, Quinn DW. An analysis of 1-D smoothed particle hydrodynamics kernels. J Comput Phys. (1996) 126:165–80. doi: 10.1006/jcph.1996.0128

CrossRef Full Text | Google Scholar

24. Panaretos VM, Zemel Y. Statistical aspects of Wasserstein distances. Annu Rev Stat Appl. (2019) 6:405–31. doi: 10.1146/annurev-statistics-030718-104938

CrossRef Full Text | Google Scholar